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ABSTRACT 


This  report  describes  23  investigations  in  the  field  of  seismic 
discrimination.  These  are  grouped  as  follows:  seismic  body 

wave  studies  (3  contributions),  seismic  surface  wave  studies 
(5  contributions),  investigations  of  seismic  source  mechanisms 
(2  contributions),  studies  of  earth  heterogeneity  and  scattering 
phenomena  (5  contributions),  studies  of  earthquake  occurrence 
rates  (3  contributions),  and  outlines  of  progress  in  seismic  data 
management  system  development  and  improvements  in  in-house 
computer  facilities  (5  contributions). 
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SUMMARY 


This  is  the  twenty-third  Semiannual  Technical  Summary  report  describing  the  activities 
of  Lincoln  Laboratory,  M,I.T. , in  the  field  of  seismic  discrimination.  These  activities  involve 
research  into  the  fundamental  seismological  problems  associated  with  the  detection,  location, 
and  identification  of  earthquakes  and  nuclear  explosions.  We  also  are  concerned  with  the  devel- 
opment of  methods  for  the  handling  and  analysis  of  large  quantities  of  global  seismic  data,  and 
the  application  of  these  methods  to  data  management  system  design  and  the  optimum  extraction 
of  scientific  information  from  high-quality  digital  data. 

Several  investigations  in  the  general  area  of  seismic  body  waves  are  described.  We  are 
continuing  to  develop  the  application  of  maximum  entropy  cepstral  analysis  to  the  detection  of 
pP-phases  for  explosions.  Delay  times  estimated  by  this  method  agree  very  well  with  published 
shot  depths  and  surface  accelerogram  recordings.  An  attempt  has  been  made  to  verify  seismic 
reciprocity  for  short -period  earthquake  data,  as  a background  to  the  problem  of  calibration  of 
aseismic  areas  in  a threshold  test  treaty  situation.  Results  are  encouraging,  but  a residual 
amplitude  variation  by  a factor  of  2 (0.3  m^  units)  remains  unexplained.  In  another  study,  a new 
improved  method  for  the  calculation  of  ellipticity  corrections  is  described. 

In  the  area  of  surface  wave  studies,  we  have  extended  previous  work  on  surface  wave  prop- 
agation in  Eurasia.  We  describe  the  first  results  from  an  extensive  attempt  to  use  existing 
crustal  structure  information  to  predict  surface  wave  path  refraction  and  multipathing  effects. 
Some  large  deviations  are  predicted.  In  another  study,  the  source  mechanisms  of  some  Hindu 
Kush  earthquakes  have  been  estimated  from  the  joint  analysis  of  spectra  from  earthquake  pairs. 
An  earth  flattening  correction  for  Rayleigh -wave  dispersion  at  long  periods  has  been  developed 
from  normal  mode  theory,  and  seems  applicable  to  a wide  variety  of  structures.  Another  inter- 
esting development  from  free  oscillation  theory  is  the  possibility  of  representing  surface  wave 
sources  in  a moment  tensor  formalism.  This  would  appear  to  be  advantageous  for  seismic  dis- 
crimination problems.  A new  approach  to  the  estimation  of  path  corrections  for  surface  wave 
magnitudes  has  been  tested  using  data  for  the  Cannikin  explosion.  Only  limited  improvement 
over  existing  techniques  has  been  demonstrated  so  far. 

Earlier  work  on  the  determination  of  seismic  source  parameters  from  long-  and  short- 
period  P-waveforms  has  been  continued.  Considerable  success  has  been  obtained  using  long- 
period  information,  and  it  is  clear  that  there  is  considerable  source  information  in  the  short- 
period  waveforms,  though  propagation  effects  make  it  hard  to  extract  in  some  cases.  An 
interesting  application  of  source  mechanism  discrimination  by  a multi-parameter  approach  is 
described.  Earthquakes  on  the  mid-Atlantic  ridge  can  be  separated  into  dip-slip  (ridge)  and 
strike-slip  (transform)  types  with  a high  level  of  success. 

Research  into  earth  heterogeneity  and  path  effects  has  concentrated  on  two  areas  — the  in- 
version of  P-wave  travel-time  residuals  to  obtain  lateral  variations  in  structure,  and  the  further 
study  of  the  scattering  of  seismic  waves.  The  first  results  from  a global  analysis  of  728,072 
observations  of  P and  PKIKP  are  presented.  Third-order  spherical  harmonic  analysis  indicates 
significant  correlations  between  inhomogeneities  at  different  levels  in  the  mantle.  Inversion  of 
travel-time  observations  at  NORSAR  shows  evidence  of  pipe-like  structures  that  extend  below 
the  Mohorovicic  discontinuity.  Polarization  studies  of  the  short -period  arrivals  at  the  three - 
component  site  at  LASA  show  a surprisingly  complex  coda  which  is  a strong  function  of  source 
location.  The  locations  of  the  scatterers  are  being  investigated.  A theoretical  study  of  primary 
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scattering  by  a transition  zone  of  changing  velocity  is  outlined.  New  attempts  to  characterize 
the  scattering  processes  beneath  large  arrays  are  described. 

We  are  continuing  investigations  into  the  properties  of  earthquake  time  series,  formed  by 
counts  of  earthquakes  above  a given  magnitude  threshold  per  unit  time.  Significant  correlations 
between  the  time  series  for  widely  separated  regions  of  various  tectonic  types  are  found.  This 
strongly  suggests  the  presence  of  a global  earthquake  triggering  mechanism.  Spectral  analysis 
has  been  applied  to  the  regional  earthquake  time  series,  using  the  maximum  entropy  method,  and 
approximately  30  percent  of  the  energy  in  the  spectrum  are  found  to  correlate  across  all  regions. 
A spectral  peak  at  a 6-month  period  is  presumably  due  to  earth  tides.  Another  at  7.6  months 
remains  unexplained. 

We  are  deeply  involved  in  the  design  and  development  of  the  data  management  system  for  the 
proposed  ARPA  seismic  network.  In  particular,  we  have  been  experimenting  with  a documenta- 
tion system  for  the  network,  based  on  the  On-Line  System  (NLS)  developed  by  Stanford  Research 
Institute.  Creation  of  file  structures  and  input  of  documentation  are  progressing.  An  attempt 
to  access  filed  information  on  the  datacomputer  is  described. 

We  continue  to  update  and  modify  our  in-house  computer  systems  so  that  we  can  take  full 
advantage  of  the  new  digital  seismic  data  as  they  become  available  on  the  datacomputer.  Some 
recent  improvements  and  progress  are  reported. 

M.  A.  Chinncry 
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GLOSSARY 


ARPA 

ARPANET 

Defense  Advanced  Research  Projects  Agency 
ARPA  Computer  Network 

ASG 

Applied  Seismology  Group 

BBN 

Bolt  Beranek  & Newman 

BMEM 

Burg  Maximum  Entropy  Method 

CCA 

Computer  Corporation  of  America 

CPU 

Central  Processing  Unit 

CRT 

Cathode  Ray  Tube 

DADS 

Data  Analysis  and  Display  System 

DDC 

Defense  Documentation  Center 

DU 

Data  Unit 

HGLP 

High  Gain  Long  Period 

ISC 

International  Seismological  Center 

IWWSS 

Integrated  World-Wide  Seismic  System 

I.ASA 

Large  Aperture  Seismic  Array 

NEIS 

National  Earthquake  Information  Service 

NLS 

On-Line  System 

NORSAR 

Norwegian  Seismic  Array 

PDE 

Preliminary  Determination  of  Epicenters 

PEM 

Parametric  Earth  Model 

SATS 

Semiannual  Technical  Summary 

SRO 

Seismic  Research  Observatory 

TENEX 

PDP-10  Operating  System 

USGS 

United  States  Geological  Survey 

WWSSN 

World-Wide  Standard  Seismograph  Network 
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SEISMIC  DISCRIMINATION 


1.  BODY  WAVE  STUDIES 


A.  MAXIMUM  ENTROPY  CEPSTRAL  ANALYSIS 


The  Maximum  Entropy  Cepstral  Analysis  method  has  been  evaluated  by  determining  the 

delay  time  between  the  P-  and  pP-phases  in  short-period  seismograms  recorded  at  Norway  for 

seven  Nevada  Test  Site  explosions  and  comparing  the  results  to  the  observed  times  determined 

2 

from  ground  zero  accelerograms.  Table  1-1  is  a list  of  the  events  and  includes  the  observed 
delay  time  and  the  calculated  times  of  the  primary  cepstral  peak  (T^,^  and  first  harmonic 
divided  by  2 {when  it  was  observed).  The  data  were  recorded  at  10  Ilz  and  the  cepstra 

were  read  to  the  nearest  0.05  sec. 


TABLE  1-1 

CALCULATED  AND  OBSERVED  DELAY 

TIMES  FOR  NTS 

EVENTS 

Date 

Origin  Time 

Event  Name 

Depth 

(m) 

TpP_p  (observed) 
(sec) 

2 

tcalc  ^calc* 

(sec) 

7-16-69 

1455:00 

Hutch 

549 

0.59 

0.45  (0.50) 

6-9-68 

1400:00 

Noggin 

582 

0.75 

0.55 

3-22-68 

1500:00 

Stinger 

668 

0.85 

0.90  (0.85) 

10-18-67 

1430:00 

Lanphere 

715 

0.86 

0.75  (0.80) 

1-19-68 

1815:00 

Faul  tless 

975 

0.80 

0.70  (0.70) 

4-26-68 

1500:00 

Boxcar 

1158 

0.96 

0.95  (0.95) 

3-26-70 

1900:00 

Handley 

1207 

0.95 

1.10 

Figure  1-1  plots  the  results.  The  observed  and  calculated  times  agree  to  within  one  or  two 
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digitizing  intervals.  Peaks  due  to  spall  closure  were  not  observed  or  possibly  were  unresolvable 
from  the  peak  due  to  the  first  harmonic  of  the  depth  phase  echo.  Since  the  existence  of  a cepstral 
peak  due  to  an  echo  depends  on  the  echo  and  primary  having  essentially  the  same  shape,  it  is  not 
surprising  that  the  peak  due  to  the  surface-spall  closure  phase  is  not  present. 

To  further  evaluate  the  method,  the  short-period  LASA  recording  of  a presumed  explosion 
from  Eastern  Kazakh  (m^  = 5.5,  A = 83°)  was  analyzed.  Figure  1-2  shows  the  seismogram  and 
the  maximum  entropy  cepstrum  indicating  a pP-P  delay  time  of  0.41  sec.  Again,  no  indication 
of  an  echo  due  to  spalling  is  observed. 

The  results  presented  here  indicate  that  the  technique  gives  accurate  estimates  of  the  pP-P 
delay  time  for  explosions  based  on  teleseismic  recordings.  Work  is  continuing  to  determine  the 
magnitude  threshold  and  minimum  delay  time  for  which  reliable  estimates  can  be  obtained. 

T.  E.  Landers 
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B.  TEST  OF  SHORT-PERIOD  RECIPROCITY 


Will  a force  F applied  at  point  P in  the  real  earth,  which  produces  a short-period  transient 
displacement  U at  point  Q,  produce  the  same  displacement  U at  P when  applied  at  Q?  This 
question  bears  on  the  feasibility  of  calibrating  certain  short-period  paths  in  a threshold  test 
treaty.  We  have  performed  an  experiment  which  attempts  to  answer  the  question  stated  above 
for  a specific  path. 

We  searched  a ten-year  earthquake  list  to  find  two  large  shallow  earthquakes  near  WWSSN 
stations  within  30°  to  90°  of  each  other.  It  developed  that  only  one  set  of  station -event  pairs  was 
found  suitable  for  the  experiment,  namely: 


Event 

20  Jan  1970 

04h  45m  00.2s 

8.8  °N  7 9.2  °W 

33  km  deep,  5.6  m^,  5.6  1V1  ^ 

4 Dec  1969 

17h  08m  48.7s 

23.1  °S  7 0.1  °W 

36  km  deep,  5.9  m^,  6.3  IV! 


Near 


Balboa  Heights,  Panama 
BHP,  8.97 °N  79.56 °W 


Antofagasta,  Chile 
ANT 


Recorded  at 


Antofagasta,  Chile 
ANT,  23.71  °S  70.43  °W 


Balboa  Heights,  Panama 
BHP 


Using  long-period  P-waves,  and  a procedure  discussed  in  Sec.  III-A  when  applied  to  shallow 
Asian  earthquakes,  source  parameters  for  each  of  these  earthquakes  were  estimated.  This  is 
not  a unique  process,  but  the  source  parameters  gained  are  consistent  with  the  long-period  data 
observed  in  the  distance  range  30°  to  90°.  We  then  used  these  source  parameters  to  generate  a 
short-period,  north-south  seismogram  at  BHP  from  the  Chilean  earthquake,  and  the  same  type 
of  record  at  ANT  from  the  Panamanian  earthquake.  We  could  not  compare  the  vertical  compo- 
nents in  this  case  since  this  record  at  ANT  from  the  Panamanian  earthquake  is  uninterpretablc. 
Because  of  the  north -south  alignment  of  the  path  considered,  the  north -south  components  are 
only  15°  off  the  radial  vector  between  the  two  sites. 

Figure  1-3  shows  the  comparison  between  the  observed  and  synthetic  records  at  BHP  and 
ANT  for  the  events  in  Chile  and  Panama,  respectively.  Although  different  source  parameters 
were  used  to  generate  the  synthetics,  the  same  spreading  factor  and  attenuation  operator  were 
applied  to  both.  The  general  form  of  the  synthetics  is  similar  to  the  observed  data  in  both  cases; 
however,  the  amplitude  of  the  synthetic  at  ANT  is  about  one-half  that  of  the  observed.  Since  the 
means  employed  here  are  somewhat  indirect,  estimating  source  parameters  in  order  to  scale 
the  amplitudes,  it  is  difficult  to  say  whether  short-period  reciprocity  failed  by  a factor  of  two 
or  our  ability  to  test  the  question  by  this  means  is  no  more  accurate  than  a factor  of  two.  The 
similarity  of  the  form  of  the  observed  and  synthetic  seismograms  in  Fig.  1-3  is,  however, 
encouraging.  j.  R Filson 

R.  G.  North 


C.  NEW  METHOD  OF  CALCULATION  OF  ELLIPTICITY  CORRECTIONS 
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The  current  usage  of  ellipticity  corrections  requires  either  extensive  tables  in  which  the 

corrections  are  listed  as  functions  of  colatitude  of  the  source,  epicentral  distance,  and  the  azi- 

4 

muth,  or  an  approximate  formula  that  may  lead  to  errors  as  large  as  0.3  sec  for  the  P-phase. 

The  simple  transformation,  invoking  the  addition  theorem  for  Legendre  polynomials,  of 
5 

Eq.  (10.30)  of  Bullen  leads  to  separation  of  variables: 


2 


(1-1) 


2 

(|  - COS2  9')  = P2°(cos  9')  = Yj  P2  m (COS  Oo)  P2  m (cose)  cos 

m=0 


where  0q  is  the  colatitude  of  the  source,  and  | is  the  azimuth  from  the  source  to  the  epicenter; 

P~>  _ are  defined: 

Z,m 


? --[e  P?m 

Z,m  m (Z  + m)IJ  Z 


c = i , e = Z for  m > 0. 
o * m 


Thus,  the  correction  can  be  sought  in  a form 
Z 

6t  = Y Tm(A)  ‘ p2  m *cos  9o'  cos  * d-2) 

m = 0 


The  efficiency  of  this  scheme  lies  in  the  fact  that  only  three  integrals  such  as  Eq.  (10.30)  of 
Bullen5  need  to  be  evaluated  for  a given  epicentral  distance  and  species  of  ray. 

We  have  investigated  the  dependence  of  ^(A)  on  the  adopted  earth  model;  two  models  were 
tested:  1066A(Ref.  6)  and  PEM  (parametric  earth  model)7  Model  1 066A  is  characterized  by 
continuous  variation  of  the  velocities  and  density  in  the  upper  mantle,  while  model  PEM  has  two 
first-order  discontinuities  at  depths  of  4Z0  and  67  0 km.  Results  for  rays  penetrating  deeper 
than  67  0 km  show  that  the  difference  in  r ( A)  for  the  two  models  is  not  more  than  0.01  sec  for 
P.  Since  the  two  models  are  thought  to  represent  two  extremes  in  modeling  the  upper  mantle, 
we  conclude  that  it  is  very  unlikely  that  any  other  acceptable  model  will  yield  significantly  dif- 
ferent values  of  Tm(A). 

Figure  1-4  shows  plots  of  r^(A)  for  the  P-phase  computed  for  surface  focus  and  a source 
at  650  km  depth.  It  has  been  commonly  assumed  that  the  effect  of  the  focal  depth  on  ellipticity 
correction  is  negligible;  yet,  in  most  unfavorable  circumstances  it  could  amount  to  0.3  sec  for 
the  deep  focus  earthquakes. 

We  also  examine  the  errors  that  could  be  introduced  by  application  of  the  approximate 
4 

formula  of  Bullen, 

6t  = (H  + h)  . f(A) 


Substituting 


H = — | e aP^>  (cos  0 ) 
3 o Z,o'  o' 

and 


2 

h = — \ e a x P~  (cos0  ) * P0  (cos  A)  • cosm£ 

3 o Aj  Z,m  o Z,m  s 

m = 0 

we  find  equivalent  functions  r ( A).  In  Fig.  1-5,  we  compare  Tm(A)  obtained  using  the  table  of 
f(A)  of  Jeffreys  and  Bullen^  with  our  computations.  The  errors  may  amount  easily  to  0.3  sec. 
The  explicit  formula  for  6t  is: 


1 Z 

6t  = ^tq(A)  (3  cos  0q  — 1 ) + *s/~3  t ^ (A ) sin  0^  cos  0^  cos  £ 


+ n/T t2(A)  sin2  0o(cos2  4 — - ) 
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All  the  trigonometric  functions  in  this  formula  are  automatically  computed  in  the  earthquake 
location  routine. 

Q 

Dziewonski  and  Gilbert  give  tables  of  rm(A)  for  eight  seismic  phases  (P,  PcP.  PKP1, 
PKP2,  PKIKP,  S,  ScS,  and  SKS)  and  three  focal  depths:  surface,  300  km,  and  650  km. 

A.M.  Dziewonski 
F.  Gilbert t 
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CEPSTRUM 


Fig.I-1.  Observed  vs  calculated  pP-P  delay  times. 
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Fig.  1-2.  Short-period  seismogram  and  maximum  entropy  complex 
cepstrum  for  presumed  Eastern  Kazakh  explosion. 
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Fig.  1-3.  Observed  and  synthetic  short-period,  north-south  seismograms  at  BHP 
(Panama)  from  earthquake  in  Chile  (left)  and  at  ANT  (Chile)  from  earthquake  in 
Panama  (right). 
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ELLIPTICITY  CORRECTION  COEFFICIENT  (sec) 


Fig.  1-4.  Comparison  of  P-phase 
ellipticity  correction  coefficients 
computed  for  focal  depths:  sur- 

face and  650  km. 


Fig.  1-5.  Comparison  of  P-phase 
ellipticity  corrections  obtained 
using  our  method  of  calculation 
t (A)  and  equivalent  values  for 
Bullen's^  approximation:  surface 
focus. 
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II.  SURFACE  WAVE  STUDIES 


A.  SURFACE  WAVE  RAY  TRACING  ACROSS  EURASIA 

Studies  of  Rayleigh-wave  multipath  propagation  using  ray  tracing  have  been  reported  for 

Western  North  America  and  the  North  America  — Pacific  boundary  by  B.  R.  Julian  in  previous 

1 

SATS.  This  work  is  now  being  extended  to  include  the  effects  of  lateral  variations  in  structure, 
particularly  those  of  major  mountain  ranges  such  as  the  Alps  and  Himalayas,  upon  the  transmis- 
sion of  surface  waves  across  Eurasia. 

Information  on  the  crustal  structure  of  this  region  is  considerably  sparser  than  for  North 
America.  With  the  exceptions  of  Western  Europe,  European  Russia,  and  the  Soviet  Far  East, 
not  much  seismic  refraction  work  has  been  carried  out  across  this  vast  area.  We  therefore  de- 
cided to  consider  initially  only  the  effect  of  variation  of  crustal  thickness  upon  the  Rayleigh-wave 

dispersion.  The  primary  source  of  this  information,  a contour  map  of  Moho  depth  in  Eastern 

2 

Europe,  the  Soviet  Union,  and  China,  is  based  (outside  the  regions  for  which  refraction  data  are 
available)  mainly  upon  elevation  with  occasional  constraints  from  gravity  observations.  This 
was  supplemented  by  similar  contour  maps  for  Scandinavia,  Western  Europe,  and  India.  No  data 
were  available  for  Spain,  Pakistan,  Afghanistan,  Iran,  the  Middle  East,  and  Southeast  Asia. 

The  area  used,  however,  constitutes  some  75  to  80  percent  of  the  Eurasian  land  mass  as  well 
as  most  of  the  adjacent  continental  shelf.  These  contour  maps  have  been  digitized  at  0.5°  inter- 
vals of  latitude  and  longitude,  and  the  resulting  grid  converted  to  one  of  phase  velocity. 

The  earth  model  chosen  for  the  dispersion  consists  of  a crustal  structure  determined  for 

3 

the  Alps  from  dispersion  data,  and  an  upper-mantle  model  calculated  to  satisfy  observed  normal 

4 

mode  eigenperiods  and  average  continental  surface  wave  dispersion.  The  effect  of  crustal  thick- 
ness has  been  approximated  by  linearly  expanding  the  crustal  layers  and  determining  the  disper- 
sion of  the  resultant  model.  Figure  II- 1 shows  the  phase  velocity  at  20-  and  40 -sec  periods  for 
Rayleigh  waves  as  a function  of  crustal  thickness.  The  values  of  phase  velocity  at  these  periods 
for  a 35-km-thick  crust  are  close  to  those  observed  for  some  Eurasian  paths.  The  phase  veloc- 
ities can  be  seen  to  increase  to  nearly  oceanic  values  as  the  crust  becomes  thinner  and  to 
decrease  quite  substantially  as  the  thickness  approaches  the  values  of  50  to  75  km  estimated  for 
the  Alps  and  Himalayas. 

Considerable  azimuthal  anomalies  have  been  observed  at  NORSAR  for  20-  and  40-sec  period 
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Rayleigh  waves,  and,  in  an  attempt  to  explain  those  observed  for  Eurasian  events,  rays  have 
been  traced  through  the  grid  at  5°  increments  from  45°  to  225°  azimuth.  Figures  II -2  and  II- 3 
show  the  results  of  this  ray  tracing  at  a 20-sec  period  on  an  azimuthal  equidistant  projection 
centered  at  NORSAR:  this  projection  has  the  advantage  that  all  great  circle  paths  through  the 
center  are  straight  lines.  The  rays  terminate  at  the  edge  of  the  grid.  Unfortunately,  because 
of  a lack  of  data  for  some  parts  of  the  Arctic  Ocean,  rays  cannot  be  traced  to  the  Kuriles  and 
Japan.  These  figures  show  clearly  the  effect  of  the  Alpine,  Himalayan,  and  Central  Asian  fold 
belts  upon  surface  wave  propagation,  as  well  as  that  of  crustal  thinning  in  continental  shelf  re- 
gions. The  number  at  the  end  of  each  ray  denotes  the  final  azimuthal  deviation  suffered. 

Figure  II-2  shows  ray  paths  traced  at  5°  from  45°  to  12  0°  azimuth  from  NORSAR.  It  can 
be  seen  that  rays  at  azimuths  up  to  75°  experience  positive  deviations;  the  ray  at  80°  is  remark- 
ably straight.  In  the  Himalayas  and  Tibet  (85°  to  105°  azimuth),  the  deviations  vary  consider- 
ably, leading  to  multipathing.  These  results  may  be  compared  with  those  measured  using  high- 
resolution  frequency-wavenumber  analysis  of  surface  wave  arrivals  at  NORSAR/*  Positive 
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(clockwise)  deviations  of  rays  traced  from  NORSAR  imply  negative  azimuthal  anomalies  for  rays 
from  an  event  to  NORSAR  and  vice  versa. 

An  event  in  Taiwan  (distance  59°,  azimuth  81°)  had  arrivals  with  azimuth  anomalies  of  about 
+10°:  this  is  unfortunately  opposite  to  those  shown.  An  event  in  Central  Russia  in  the  same  di- 
rection gave  a similar  anomaly,  but  later  arrivals  had  negative  values  which  are  more  in  agree- 
ment with  those  predicted.  Early  arrivals  from  events  in  the  Hindu  Kush  (45°,  90°)  and  in  Tibet 
(59°,  81°)  came  in  with  small  negative  anomalies;  later  ones  were  highly  positive.  The  model 
predicts  both  these  effects,  though  the  deviations  are  about  a factor  of  2 to  3 times  less  than  the 
observed  anomalies.  The  velocity  model  evidently  needs  considerable  refinement  to  produce  the 
very  high  values  observed. 

Figure  II- 3 shows  rays  traced  at  5°  increments  from  12  5°  to  22  5°  azimuth  from  NORSAR. 
Although  crustal  thicknesses  in  the  Alps  approach  those  for  the  Himalayas,  the  effect  is  not  as 
marked,  primarily  because  the  Alps  run  normal  to  the  ray  paths  in  most  cases.  The  effect  of 
the  Carpathians  is  somewhat  greater,  and  it  can  be  seen  that  there  is  considerable  divergence 
of  the  ray  paths  in  the  Balkans.  This  should  result  in  higher-than-average  surface  wave  ampli- 
tudes at  NORSAR  for  events  in  the  Aegean  and  Western  Turkey:  a hypothesis  which  we  hope  to 
verify  by  observation.  Fairly  large  deviations  also  occur  around  the  coast  of  Denmark.  Mea- 
surements of  arrival  azimuth  for  surface  waves  approaching  NORSAR  from  the  south  show  fairly 
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large  deviations,  attributed  to  refraction  at  the  coastline  of  North  Africa.  The  present  results 
indicate  that  at  least  part  of  this  (up  to  15°  deviation)  may  be  caused  by  refraction  at  structures 
within  20°  distance  from  NORSAR. 

This  study  is  going  to  be  extended  to  include  variations  in  crustal  structure,  and  we  hope 
that  this  will  help  to  explain  the  much  larger  deviations  observed.  A division  into  shields,  plat- 
forms, and  foldbelts  should  give  greater  variations  in  phase  velocities.  The  Caspian  and  Black 
Seas  are  believed  to  contain  fragments  of  oceanic  crust,  and  inclusion  of  such  features  should 
produce  similar  variations.  An  extension  of  the  area  studied  to  oceanic  regions  may  be  fairly 
easy  since  there  seem  to  be  simple  relations  between  age  zones  and  spreading  rates,  and  litho- 
spheric thickness,  sedimentary  cover,  and  water  depth.  By  summing  the  group  delay  along  the 
ray  path,  it  also  should  be  possible  to  predict  the  arrival  time  of  each  refracted  wave. 

R.  G.  North 

B.  AN  EMPIRICAL  EARTH  FLATTENING  CORRECTION 

FOR  RAYLEIGH-WAVE  DISPERSION 

Calculations  of  dispersion  curves  for  a layered  half-space  using  Haskell*  s^  approach  are 

simpler,  faster,  and  require  less  computer  memory  than  computations  for  the  spherical,  grav- 
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itating  earth.  For  the  Love  waves,  Gerver  and  Kazhdan  and  Biswas  and  Knopoff  proposed  an 
exact  earth  flattening  transformation.  A similar  transformation  for  Rayleigh  waves  is  much 
more  difficult,  if  not  impossible,  to  achieve.  Bolt  and  Dorman  derived  an  empirical  formula 
for  correcting  the  Ravleigh-wave  phase  velocities  obtained  for  a flat  earth  model,  and  their  for- 
mula commonly  has  been  used  for  the  last  10  to  15  years.  Because  of  the  progress  in  the  pre- 

10 

cision  and  stability  of  the  numerical  solutions  to  the  normal  mode  problem,  we  decided  to  re- 
examine the  question  of  the  earth  flattening  correction. 

4 

We  have  used  an  oceanic  (PEM-O)  and  a continental  (PEM-C)  model  of  the  earth.  These 
models  satisfy,  in  the  sense  of  a weighted  average,  the  observed  periods  of  free  oscillations  and 
typical  surface  wave  dispersion  for  each  type  of  region.  Both  models  are  identical  below  420  km 
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depth.  As  the  flat  earth  program  requires  a discrete  layered  representation  of  the  earth,  we 
have  introduced  appropriate  modifications  in  the  PEM  models. 

The  spherical  earth  normal  mode  eigenfrequencies  are  computed  with  and  without  con- 
sideration of  gravity  for  the  modes  (Rayleigh  wave).  The  phase  velocity  for  a spherical 
earth  is 

Cs  = a • oU|/(l  + 1/2) 

a being  the  radius  of  the  earth.  The  phase  velocities  in  a flat  layered  half-space  C,  have  been 

£ I 

computed  at  frequencies  using  a program  based  on  a layer  matrix  formulation.  Results  ob- 
tained using  this  program  agree  to  the  published  accuracy  (usually  5 to  6 significant  figures) 

1 1 

with  those  of  other  matrix  algorithms  and  also  with  those  of  a program  based  on  a numerical 
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integration  method. 

Consideration  of  gravity  in  the  normal  mode  solutions  has  the  effect  of  increasing  the  phase 
velocity.  This  increase  is  0.0002  km/sec  at  the  50-sec  period,  rising  to  0.0032  km/sec  at 
150  sec  and  to  0.021  at  the  300-sec  period.  As  the  effect  of  gravity  is  only  0.005  percent  at 
50  sec,  calculations  for  shorter  periods  are  carried  out  for  a nongravitating  earth.  In  general, 
the  effect  of  gravity  in  the  period  range  considered  is  an  order  of  magnitude  less  than  that  of 
sphericity:  at  300  sec  6 C = Cg  — Cf  is  0.24  km/sec;  for  the  continental  model  6 C is  still  as  large 
as  0.01  km/sec  at  the  25-sec  period,  and  for  the  oceanic  model  this  difference  extends  to  periods 
as  short  as  15  sec. 

In  Fig.  II-4,  a plot  of  6C/T  (division  by  period  for  scaling  purposes)  shows  a variation  which 
in  shape  resembles  much  more  closely  the  group  velocity  variation  (chain  line)  than  the  phase 
velocity  curve  (broken  line).  Like  the  group  velocity  curves  <5C/T  possesses  a short-period 
maximum  and  a long-period  minimum,  although  the  minima  do  not  coincide  precisely. 

9 

The  line  0.00016  • is  equivalent  to  the  Bolt-Dorman  correction.  Having  noted  the  sim- 
ilarity between  6C/T  and  u^,  we  determine  the  proportionality  factor: 

n = (C  - C-)/(u.  • T) 
u ' s f ' ' f 

For  the  continental  model,  the  best  fitting  value  of  qu  is  0.000156,  and  for  the  oceanic  structure 
it  is  0.000166;  the  mean  of  the  two  (0.000162)  still  provides  a reasonably  good  fit. 

In  Fig.II-5,  we  compare  errors  associated  with  our  empirical  correction  (<5C  = a • uf  • T) 
with  that  of  Bolt  and  Dorman.  For  both  formulas,  the  errors  are  very  small  for  periods  less 
than  100  sec;  this  means  that  systematic  bias  introduced  by  application  of  either  of  the  two  em- 
pirical corrections  is  insignificant  in  comparison  with  measurement  errors.  Above  100  sec, 
however,  the  advantage  of  using  6C  = <*u  • u^  • T is  distinct;  it  predicts  the  exact  values  of  Cg 
to  within  0.2  percent  (less  than  0.01  km/sec)  to  periods  as  long  as  230  sec,  which  is  the  approx- 
imate upper  limit  of  the  minor  arc  measurements  of  surface  wave  dispersion. 

Since  the  same  value  of  fits  the  two  quite  different  models  fairly  well,  it  seems  that  the 
above  relation  will  be  a good  approximation  for  a wide  variety  of  structures. 

R.  G.  North 
A.M.  Dziewonski 
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C.  MOMENT  TENSOR  REPRESENTATION  OF  SURFACE  WAVE  SOURCES 


The  usual  equations  for  Love  and  Rayleigh  waves  due  to  a point  source  in  a layered  elastic 
12 

half-space  (e.g.,  Saito  ) can  be  expressed  as  a linear  combination  of  the  moment  tensor  elements. 

The  advantages  of  such  an  approach  are  (1)  that  the  expansion  is  linear  and  therefore  amenable 
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to  the  stochastic  least-squares  estimation  procedure  (Foster  ),  and  (2)  that  the  expansion  is 

14 

complete  consisting  of  monopole,  dipole,  and  quadrupole  terms  (Knopoff  and  Randall  ) which 

ensures  reliability  in  the  estimated  source  parameters.  An  easy  way  to  prove  this  expansion  is 

15 

to  merge  the  moment  tensor  formalism  of  Gilbert  with  the  vector  harmonic  expressions  given 
by  Saito.*2 

The  Love  and  Rayleigh  wave  problems  separate  out  onto  vector  components  when  written  in 
terms  of  the  vector  cylindrical  harmonics.  These  are: 
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where  is  the  outward  traveling  Hankel  wave  function.  The  expansion  becomes: 
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where  the  expansion  coefficients  1,  r^,  and  r^  are  solutions  to  the  well-known  Love  and  Rayleigh 
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Introducing  expansion  weights  a^  and  b^  allows  the  expansion  to  be  written: 
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where  the  sum  over  j is  over  the  individual  Love  and  Rayleigh  wave  modes  of  propagation. 

Since  w only  appears  in  the  expansion  weights,  the  expression  can  be  integrated  for  a step 
function  source  time  dependence  as  Gilbert*5  did  for  the  free  oscillation  expansion  to  give  exact 
expressions  for  aj  and  b^.  These  are: 
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where  the  normalization  integral  is 
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where  the  normalization  integral  is 

r00  2 2 

IR(w.j<k)5j  dzp(z)  [|r1(«j,k,  z)|  + | k,  z)|  ] 

In  both  Eqs.(II-S)  and  (II-6),  <5  is  the  Kronecker  delta  and  zg  refers  to  the  source  depth. 

To  recapitulate,  Eqs.  (II- 5 ) and  ( II- 6 ) imply  the  following  assumptions: 

(1)  The  source  and  receiver  exist  in  the  same  layered  half-space, 

(2)  The  source  exists  at  only  one  point  in  the  half-space, 

(3)  The  source  has  a step  function  time  dependence, 

(4)  All  body  wave  effects  are  ignored,  and 

(5)  The  material  is  unstressed  before  the  source  acts. 

With  these  assumptions,  the  six  elements  of  the  moment  tensor  provide  a complete  description 
of  the  point  source.  It  can  be  decomposed  uniquely  into  three  source  mechanisms:  (1)  symmet- 
ric or  explosive  source,  (2)  compensated  linear  vector  dipole  (CLVD)  source,  and  (3)  double 

couple  source,  provided  that  the  principal  stress  directions  for  the  latter  two  are  coincident 
14 

(Knopoff  and  Randall  ).  As  such,  the  expansion  is  a good  candidate  for  the  application  of  linear 
estimation  procedures  and  should  provide  discriminant  information  between  explosive  and  earth- 
quake sources.  _ ,T.  ~ 

M D.  W.  Me  Cow  an 
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D.  TEST  OF  THE  EMPIRICAL  SURFACE  WAVE  MAGNITUDE  PATH  CORRECTION 

1 6 

In  the  last  SATS,  an  empirical  method  of  computing  the  path  correction  for  surface  wave 
magnitudes  (Mg)  was  suggested.  This  method  is  based  on  the  saddle-point  approximation  for 
the  amplitudes  of  damped  surface  waves.  One  term  of  this  approximation  is  inversely  related 
to  the  slope  of  the  group  velocity  curve,  a term  that  had  heretofore  been  ignored  (or  taken  from 
tables  for  specific  paths)  in  the  computation  of  Mg.  The  crux  of  the  proposed  method  is  to  esti- 
mate the  slope  of  the  group  velocity  curve  directly  from  the  individual  seismogram  on  which  the 
surface  wave  amplitudes  are  measured.  This  estimate  entails  measuring  the  period  and  arrival 
time  of  the  surface  wave  phases  arriving  just  before  and  after  the  phase  on  which  the  amplitude 

to  be  used  in  the  M calculation  is  measured, 
s 

This  method  has  been  tested  using  the  long-period  vertical  seismograms  recorded  at  some 
thirty  World-Wide  Standard  Seismograph  Network  (WWSSN)  stations  following  the  large  under- 
ground nuclear  explosion  "Cannikin,"  detonated  in  the  Aleutian  Islands  on  6 November  1971.  We 
decided  to  test  the  method  with  an  explosive  source  in  order  to  minimize  the  azimuthal  variation 
in  surface  wave  amplitudes  predicted  from  the  double  couple  earthquake  source.  The  "new"  Mg 
formula  is  based  on  the  logarithm  of  the  amplitude  of  the  surface  wave  phase,  plus  correction 
terms  for  geometric  spreading,  attenuation,  and  the  group  velocity  slope  estimation.  The  re- 
sults were  compared  with  the  "old"  1V1  g (or  Prague)  formula  which  is  based  on  the  sum  of  the 
logarithm  of  the  ratio  of  the  amplitude  and  period,  and  a geometric  spreading  term.  A positive 
result  from  the  experiment  would  be  that  the  new  method  would  reduce  significantly  the  scatter 
(or  standard  deviation)  of  Mg  measurements,  at  various  periods  and  distances,  over  the  old 
method. 

In  Fig.  II-6,  the  old  and  new  methods  of  Mg  determinations,  as  a function  of  period,  are  com- 
pared using  the  Cannikin  data.  In  this  and  the  following  comparisons,  only  the  variable  factors 
affecting  Mg  are  considered;  the  addition  of  an  arbitrary  constant  (3.3  units  in  the  case  of  the 
old  measurement)  has  not  been  made.  Figure  II- 6 shows  that,  although  the  new  Mg  does  not  show 
the  dependence  on  period  seen  in  the  old,  the  overall  scatter  in  the  measurements  does  not  seem 
to  have  been  significantly  reduced.  The  trend  in  the  old  measurements  can  be  attributed  almost 
wholly  to  the  use  of  the  ratio  of  amplitude  and  period  rather  than  amplitude  alone. 

Figure  II- 7 compares  the  old  and  new  methods  as  a function  of  distance.  In  this  diagram, 
it  is  generally  the  case  that  at  a single  station  (distance  = constant)  the  new  measurements  are 
less  scattered  than  the  old;  however,  the  overall  scatter  is  about  the  same.  This  means  that  the 
new  method  will  give  more  consistent  results  using  phases  of  different  periods  at  a single  sta- 
tion, but,  again,  the  scatter  between  various  phases  at  different  stations  is  essentially  the  same. 

These  conclusions  are  quantified  in  the  histograms  of  Fig.  11-8  which  shows  the  distribution 
of  the  old  and  new  measurements.  The  standard  deviation  of  the  old  method  is  0.24  unit;  that  of 
the  new  is  0.22.  This  does  not  appear  to  be  an  improvement  of  enough  significance  to  encourage 
the  general  use  of  the  proposed  new  method  of  surface  wave  magnitude  determination.  Three 
reasons  for  this  lack  of  improvement  are  initially  suggested.  An  assumed  attenuation  correction 
based  on  a constant  Q of  500  was  applied  in  the  new  method.  Both  azimuthal  and  frequency  de- 
pendence of  Q may  make  this  assumption  invalid.  Most  paths  from  the  Aleutians  to  any  continen- 
tal station  are  complicated  by  crossing  ocean-continent  boundaries.  This  may  introduce  errors 
which  cannot  be  accounted  for  in  the  simple  propagation  model  applied.  Finally,  it  may  be  that 
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differences  in  arrival  time  and  period,  used  in  the  estimation  of  the  group  velocity  curve  slope, 
cannot  be  measured  accurately  enough  by  hand. 

J.  R.  Filson 

E.  RAYLEIGH -WAVE  PHASE  VELOCITIES  ACROSS  THE  EURASIAN  CONTINENT 

Four  earthquakes  located  within  a small  region  of  the  Hindu  Kush  have  been  studied  thor- 
oughly to  determine  their  source  mechanisms  (depth,  slip  angle,  dip  angle,  strike,  and  seismic 

1 6 

moment).  Two  of  these  earthquakes  (Nos.  1 and  2)  were  analyzed  in  a previous  report  using 

1 7 

a surface  wave  method  developed  by  Weidner  and  Aki.  This  method,  which  is  restricted  to 
pairs  of  earthquakes,  recently  has  been  extended  to  incorporate  any  number  of  events  in  a clus- 
ter. Furthermore,  highly  refined  Rayleigh-wave  phase  velocities  for  paths  emanating  from  the 
source  region  to  surrounding  WWTSSN  stations  have  been  obtained. 

Unlike  Weidner  and  Aki’s  method  involving  spectral  ratios  between  an  earthquake  pair,  the 
new,  extended  method  isolates  source  effects  on  surface  wave  spectra  for  a single  earthquake 
and  then  matches  theoretical  calculations  to  obtain  source  parameters.  To  accomplish  this,  we 

use  earthquakes  with  known  source  mechanisms  to  estimate  the  propagational  effects-amplitude 

th 

attenuation,  IL(o;),  and  path  phase  delay,  (p ^(co ),  between  the  source  region  and  the  j station. 

Then,  we  deconvolve  these  effects  from  spectra  of  a new  earthquake,  thereby  isolating  its  source 

effects.  For  example,  if  E.  .(a;)  is  the  known  amplitude  excitation  for  earthquake  No.  1 and  A.  .(t o) 

th 

is  the  observed  amplitude  spectrum  at  the  j station,  corrected  for  instrument  response  and 
geometric  spreading,  then 

Ay(0>  ) 


J 

Similarly,  another  estimate  for  II. (a; ) can  be  found  using  E-.(u;)  and  A-.{co).  Since  these  estimates 

J J ^ g J 

are  subject  to  errors,  we  make  use  of  a statistical  calculation  to  obtain  the  maximum-likelihood 
estimate  (MLE),  IL(u;).  Deconvolving  this  estimate  from  spectra  of  a new,  third  earthquake  oc- 
curring in  the  source  region 

A (a;) 

S3  4<«>  = -f1- 

-1  H.(«) 


we  isolate  its  source  effects,  which  are  then  modeled.  It  should  be  noted  that  (1)  an  anal- 

ogous procedure  is  carried  out  on  the  phase  spectra  involving  the  MLE  for  the  path  phase  delay, 
<pj( u),  and  that  (2)  as  more  events  in  the  source  region  are  solved,  the  estimates  IL(o;)  and  <p^(oo) 
are  steadily  improved. 

Using  <pj {co)  obtained  after  four  events  were  solved,  wc  calculated  Rayleigh-wave  phase  ve- 
locities for  many  paths  between  the  source  region  and  the  stations  shown  in  Fig.  II-9.  The  dis- 
persion for  just  the  paths  which  sample  the  stable  Russian  platform  are  shown  in  Fig.  11-10.  The 

dotted  lines  in  the  first  group  indicate  the  range  of  dispersion  observed  for  stable  areas  of  the 

1 9 

North  American  continent  compiled  by  Knopoff.  The  next  group  samples  the  Norwegian  Sea  as 

well  as  the  platform,  and  thus  these  velocities  are  noticeably  slower.  Although  the  last  group 

shows  some  scatter  at  long  periods,  comparison  with  the  first  group  is  good  at  periods  shorter 

than  70  sec.  „ T 

II.  J.  Patton 
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Fig.  II— 1 . Rayleigh-wave  phase  velocities  at  20-  and  40-sec  periods 
for  crustal  model  used,  as  a function  of  crustal  thickness. 
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Fig.II-2.  Ray  paths  traced  from  NORSAR  (60.8°N,  10.8°E)  at  5°  increments 
in  azimuth  from  45°  to  120°.  Numbers  at  end  of  each  ray  give  final  azimuthal 
deviation.  Longest  ray  is  to  78°  distance  from  NORSAR. 


19 


Fig.  II— 3 . Same  as  Fig.  II-2,  but  from  12  5°  to  22  5°  azimuth.  Longest  ray  is  to  29° 
distance  from  NORSAR. 
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Fig.  II— 4 , Comparison  of  computed  values  of  (Cs  — Cf)/T  with  scaled  phase 
(Cf)  and  group  (uf)  velocity  curves  computed  for  layered  half-space. 
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Fig.  II-5.  Comparison  of  relative  errors  introduced  by  application  of  empirical 
earth  flattening  correction  proposed  in  this  study  (solid  line)  with  that  of  Bolt  and 
Dorman.^  designates  value  of  phase  velocity  obtained  using  the  appropriate 
approximate  formula;  C-g  is  the  exact  value. 
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Fig.  II-6.  Comparison  of  new  and  old  Ms  measurement  methods  as  a 
function  of  period. 


Fig.  II- 7 . Comparison  of  new  and  old  Ms  measurement  methods  as  a 
function  of  distance. 
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Fig.  II-8.  Distribution  of  new  and  old  Ms  measurement  methods. 
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Fig.  II— 9 • Station  distribution  used  in  this  study. 
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Fig.  II - 1 0 . Dispersion  curves  for  paths  that  sample  the  Russian  platform. 
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III.  EARTHQUAKE  SOURCE  MECHANISMS 


A.  P-WAVES  FROM  SHALLOW  ASIAN  EARTHQUAKES 

1 

As  indicated  in  the  last  SATS,  a study  into  the  usefulness  of  long-  and  short-period 

P-waveforms  in  determining  source  parameters  and  depth  of  shallow  focus  earthquakes  has 

been  undertaken.  To  date  some  seven  shallow,  rather  large  events  from  Asia  have  been  studied 

using  long-period  data.  The  method  of  study  consists  of  visually  matching  the  long-period 

P-waves  observed  at  a network  of  stations  (here  the  WWSSN)  to  those  generated  by  a simple 

model.  The  model  seismograms  are  based  on  the  sum  of  the  phases  P,  pP,  and  sP  generated 

2 

by  a Savage-type  source,  this  sum  then  being  convolved  with  a nominal  instrument  response, 
an  attenuation  operator,  and  a geometric  spreading  factor.  Account  is  made  for  the  amplitude 
and  phase  modification  of  pP  and  sP  due  to  reflection  above  the  source,  and  for  the  effect  of  the 
free  surface  at  the  recording  site.  The  effects  of  crustal  layering  near  the  source  and  receiver 
are  ignored. 

A sample  of  the  results  of  this  type  of  analysis  are  given  in  Fig.  Ill— 1 through  the  example 
of  a shallow  earthquake  in  Kashmir.  Here,  13  observed  seismograms  (above)  may  be  compared 
with  the  s^Tithetic  seismograms  generated  for  that  station  (directly  below)  by  the  model.  Both 
the  data  and  the  synthetic  seismograms  have  been  scaled  to  long-period,  WWSSN,  vertical  seis- 
mograms at  a gain  of  1500.  The  scale  gives  the  equivalent  time  and  amplitude  on  such  an  in- 
strument; the  model  parameters  are  shown  beneath  the  diagram.  A and  B are  the  semimajor 
and  minor  axes  of  a plane,  elliptical  fault  surface;  RVEL  is  the  velocity  of  rupture  from  one 
focus  of  the  ellipse;  C is  the  P-phase  velocity  in  the  source  region;  STRIKE  is  that  of  the  fault 
plane  measured  clockwise  from  north;  TILT  is  the  angle  the  fault  plane  makes  with  the  vertical; 
SLIP  and  DLOC  define  the  direction  and  magnitude  of  a constant  dislocation  vector  within  the 

fault  plane;  and  M gives  the  moment  based  on  the  product  of  A • B • DLOC  and  a rigidity  of 
11  2 ° 

3 • 10  dyn/cm  . The  circle  in  the  center  of  the  figure  represents  a stereographic  projection 
of  the  lower  focal  sphere  showing  the  fault  and  auxiliary  planes,  and  the  point  at  which  the  rele- 
vant rays  penetrate  this  hemisphere.  Here,  the  plane  on  the  upper  right  is  taken  to  be  the  fault 
plane.  In  general  it  may  be  said  that  the  model  gives  a fair,  visual  representation  of  the  observed 
data,  both  in  amplitude  and  phase,  for  the  first  one-and-one-half  cycles.  A similar  example  of 
another  shallow,  thrust  event,  this  time  from  the  Hindu  Kush,  is  shown  in  Fig.  Ill  — 2 . 

An  attempt  to  apply  the  same  technique  to  a strike-slip  source  is  shown  in  Fig.  Ill— 3 . This 
event,  in  northeastern  Siberia,  has  been  discussed  previously.  As  the  reader  may  observe, 
this  match  between  the  recorded  and  synthetic  seismograms  is  not  as  good  as  in  the  case  of  the 
thrust-type  sources  discussed  above.  However,  the  first  cycle  or  so,  again  in  both  phase  and 
amplitude,  is  fairly  well  reproduced.  The  source  parameters  used  to  fit  the  long-period  seis- 
mograms in  Fig.  Ill -3  were  then  used  to  generate  short -period  traces,  and  four  of  these  are  com- 
pared with  the  relevant  observations  at  four  stations  in  Fig.  Ill -4.  It  would  seem  that  this  is  the 
most  important  question  that  studies  such  as  this  address  — can  reliable  source  information  be 
gained  from  short-period  seismograms?  A comparison  of  the  observed  and  synthetic  traces  in 
Fig.  Ill -4  shows  that  the  short-period  synthetics  give  a fair  approximation  to  the  general  shape  of 
the  observed  waveforms.  Other  comparisons,  not  shown  here,  are  very  poor.  The  conclusion, 
thus  far,  is  that  along  certain  paths  information  concerning  the  earthquake  source  is  discernible 
in  the  short-period  seismogram,  while  along  other  paths  this  information  is  totally  masked  by 
propagation  effects.  j.  R Filson 
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B.  EARTHQUAKE  SOURCE  MECHANISM  DISCRIMINATION 


Source  mechanism  discrimination  for  earthquakes  of  the  mid -Atlantic  ridge  system  has 
been  studied  and  shown  to  be  feasible.  This  region  was  studied  because  the  mechanisms  of  mid- 
ocean ridge  earthquakes  are  characteristically  of  two  types:  strike -slip  events  occurring  on 

3 4 

transform  faults,  and  dip-slip  events  occurring  on  ridge  crests.  ' In  addition,  since  the  earth- 

5 

quakes  of  these  two  mechanism  types  occur  at  similar  depths,  the  study  of  the  two  mechanism 
types  can  be  approached  as  a two-class  discrimination  problem. 

Parameters  of  the  teleseismic  signal  other  than  its  first  motion  were  studied  in  order  to 
identify  those  parameters  of  the  signal  which  reflect  the  source  mechanism  differences.  Among 
the  parameters  considered  were  surface  and  body  wave  magnitudes,  P-wave  signal  complexity, 
dominant  period,  and  spectral  ratio.  A set  of  prototype  events  of  strike -slip  (SS)  or  dip-slip 
(DS)  mechanism  was  studied  to  determine  which  of  these  parameters  do  in  fact  exhibit  differences 
for  the  two  mechanism  types.  Those  parameters  found  useful  for  the  mechanism  discrimination 
were  Mg,  m^,  and  the  predominant  period  of  the  P-wave  (at  LASA).  Features,  or  functions 
having  characteristically  different  values  for  events  of  the  two  classes,  were  then  formed  from 
these  parameters,  and  the  discrimination  was  based  on  these  features. 

The  surface  and  body  wave  magnitudes  for  the  events  of  the  central  mid -Atlantic  ridge 
(from  0°  to  35°N)  studied  here  are  plotted  in  Fig.  Ill  - 5.  Similarly,  the  P-wave  predominant 
period  for  these  events  is  plotted  against  surface  wave  magnitude  in  Fig.  Ill -6.  The  two  popu- 
lations are  reasonably  separated  in  these  figures,  and  possible  lines  of  separation  have  been 
drawn.  The  features  for  the  mechanism  discrimination  are  based  on  projection  of  the  observed 
data  onto  the  normals  to  these  separation  lines.  The  performance  of  this  discrimination  scheme 
is  indicated  by  the  confusion  matrix: 


SS 

DS 

Total 

SS 

24 

1 

25 

DS 

1 

11 

12 

Thus  for  these  data,  the  scheme  yields  one  error  out  of  25  events  in  the  recognition  of  SS  events, 
and  one  error  out  of  12  events  in  the  recognition  of  DS  events.  An  estimate  of  the  total  probable 
error  rate  is  perhaps  12  percent,  or  an  88-percent  correct  classification  rate.  However,  this 
type  of  error  estimate,  based  on  the  prototype  event  set,  is  at  best  a lower  bound  to  the  true 
probable  error f*  and  must  be  further  qualified  by  the  relatively  small  sample  size.  However, 
the  separation  of  the  populations  in  Figs.  Ill- 5 and  III- 6,  on  which  the  discrimination  is  based, 
is  statistically  significant  at  the  95-percent  level. 

These  results  can  be  interpreted  in  terms  of  the  physical  environment  at  the  ridge.  The 
differences  which  might  be  observed  in  M and  m^,  due  solely  to  differences  in  the  surface  and 

body  wave  radiation  patterns  for  different  types  of  faulting,  can  be  estimated  from  the  formula- 
7 8 

tions  of  Harkrider  and  Ben-Menahem  et  aL  For  the  same  moment  and  source  time  function, 
the  estimated  changes  in  Mg  and  m^  from  their  values  for  vertical  strike -slip  faulting  are,  for 
vertical  dip-slip  faulting  and  dip-slip  faulting  on  a 45°  plane: 

AM  Am, 

s b 

Vertical  DS  —0.2  +0.6 

45°  DS  -0.3  +0.8 
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In  Fig.  Ill -5,  the  Mg  vs  positions  of  a vertical  dip-slip  event  and  a 45°  dip-slip  event  are 
indicated,  given  that  a vertical  strike-slip  event  has  an  Mg  vs  typical  for  the  observed  SS 
population. 

The  presence  of  a zone  of  high  attenuation  under  ridge  crests  but  not  under  transform 
faults  will  have  an  effect  on  the  Mg-m^  values  opposite  to  that  of  the  source  mechanisms.  The 
source  mechanism  effect  will  tend  to  separate  the  two  populations,  while  a zone  of  high  attenua- 
tion under  the  ridge  DS  events  will  reduce  their  m^,  driving  the  DS  population  back  into  the  SS 
population.  If  the  clustering  of  the  two  populations  in  Fig.  Ill -5  is  due  to  greater  attenuation  for 
DS  events,  then  that  attenuation  is  apparently  causing  a reduction  of  the  DS  body  magnitudes  by 
roughly  0.5  unit  more  than  the  reduction  in  SS  body  magnitudes.  If  the  transform  events  are 
relatively  unattenuated,  and  the  attenuating  zone  under  the  ridge  crest  is  of  the  dimensions 
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(~  1 00  km)  suggested  by  Solomon  and  Julian,  then  the  reduction  in  magnitude  requires  a Q of 

10 

about  40  in  this  zone,  consistent  with  the  estimate  by  Solomon  of  Q ~ 10  for  S-waves.  Also, 
if  the  clustering  of  the  populations  is  due  to  greater  attenuation  for  the  DS  events,  discrim- 
ination between  these  two  mechanism  types  should  be  easier  in  a region  without  this  differential 
attenuation.  A.  Shakal 
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KASHMIR  3 SEPT  1972  6.3  6.2  DEPTH  = 10km 

A = 8 km  STRIKE  * 310°  MQ=  9.0  1025 

B = 8 km  TILT  = 30° 

RVEL  = 2. 5 km/sec  SLIP  2 -90" 

C « 6.0  km/sec  DLOC  s 1.5  m 


Fig.  Ill— 1 . Comparison  of  observed  (above)  and  synthetic  (below)  long-period  P-wave 
seismograms  for  shallow  thrust  earthquake  in  Kashmir. 
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HINDU  KUSH  24  JUNE  1972  6.0  6.1 

A = 8km  STRIKE  = 25°  DEPTH  = 16km 

B * 8km  TILT  = -25°  M0  - 9.0  1025 

RVEL  * 2.5km/sec  SLIP  = 90° 

C=6.0km/sec  DLOC  * 1.5m 


Fig.  III-2.  Comparison  of  observed  (above)  and  synthetic  (below)  long-period  P-wave 
seismograms  for  shallow'  thrust  earthquake  in  Hindu  Kush. 
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SIBERIA  18  MAY  1971 
A = 16  km 
0 s 16  km 

RVEL  = 2.5  km/sec 
C = 6.0  km/sec 


5.8  6.6 

STRIKE  = 313°  DEPTH  = 12  km 
TILT  * 10°  M0  - 4.8-1026 

SLIP  = -5° 

DLOC  * 2.0m 


Fig.  Ill  - 3.  Comparison  of  observed  (above)  and  synthetic  (below)  Long-period  P-wave 
seismograms  for  shallow  strike-slip  earthquake  in  Siberia. 
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DELTA  : 64.43 
AZID  = 311.7 
TOA  : 20.7 
TSTAR  - 0.40 
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Fig.  Ill -4.  Comparison  of  observed  and  synthetic  short -period  seismograms 
at  four  European  stations  for  Siberian  event  of  Fig.  Ill -3. 
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Fig.  Ill -5.  Surface  and  body  wave 
magnitudes  for  central  mid-Atlantic 
ridge  events  of  this  study,  with  pos- 
sible line  of  separation  between  two 
event  types.  Also  indicated  are  rel- 
ative Ms-vs-m^  positions  of  90°  I)S 
and  45°  DS  events,  estimated  in  text, 
given  that  Ms  vs  mb  of  a vertical  SS 
event  is  as  indicated. 


Fig.  Ill -6.  Predominant  period  of  LASA 
P-wave  vs  Ms  for  events  of  Fig.  Ill  - 5 for 
which  LASA  signal  was  available. 
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IV.  EARTH  HETEROGENEITY  AND  SCATTERING 


A.  SYSTEMATIC  INVERSION  FOR  LARGE-SCALE  HETEROGENEITIES 

IN  THE  MANTLE-  PRELIMINARY  RESULTS 

Global  studies  of  the  radial  velocity  distribution  have  been  limited  so  far  to  the  "average" 
earth.  This  corresponds  to  the  zero-order  coefficient  in  the  spherical  harmonic  expansion  of 
the  velocities  in  the  laterally  heterogeneous  earth.  Availability  of  a large  number  of  travel 
time  observations  reported  to  the  International  Seismological  Center  (ISC)  makes  feasible  an 
attempt  at  extracting  the  coefficients  for  several  low-order  harmonics  of  the  lateral  variations 
in  the  distribution  of  compressional  velocity.  If  the  experiment  were  successful,  some  funda- 
mental geophysical  questions  might  be  answered.  Also,  availability  of  a three-dimensional 
velocity  model  could  improve  the  accuracy  of  location  of  seismic  events. 

Data  for  P and  PKIKP  travel  times  were  extracted  from  the  ISC  Bulletins  for  the  years 

1964  to  1970  for  events  with  magnitude  ^4.8  or  more  with  50  stations  reporting.  Data  for  phases 

that  are  subject  to  frequent  misinterpretations  were  rejected,  particularly  P for  distances 

<27°  and  PKIKP  for  A < 160°.  Also,  arrival  times  with  residuals  >5  sec  were  discarded  as 

1 

gross  errors.  After  correction  for  ellipticity,  the  total  number  of  retained  observations  is 

2 

728,072.  Differences  between  the  ISC  travel  times  and  the  times  computed  for  the  PEM  rep- 
resent the  data  set  used  to  determine  the  three-dimensional  perturbations  in  velocity. 

The  inverse  problem  is  approached  by  the  linear  approximation  of  the  perturbation  of  travel 

3 1 

times  by  small  changes  in  velocity  distribution.  * A discrete  model  is  adopted  in  this  study, 
although  other  representations  also  are  feasible. 

Discretization  is  achieved  by  dividing  the  mantle  into  several  concentric  shells,  and  then 
partitioning  each  shell  into  conical  blocks  of  equal  latitudinal  and  longitudinal  extent.  Pertur- 
bation of  velocity  is  assumed  to  be  constant  within  a block.  During  its  travel  through  each  con- 
centric shell,  a seismic  ray  may  cross  one  or  more  boundaries  between  adjacent  blocks;  in  this 
case,  it  is  assumed  that  the  contribution  to  the  observed  travel  time  residual  comes  from  the 
block  in  which  the  ray  spends  the  most  time.  An  example  of  actual  discretization,  together 
with  results  of  inversion,  is  shown  in  Table  IV-1. 

The  matrix  inversion  is  performed  using  the  singular  value  decomposition  technique.  For 
the  example  shown,  119  out  of  120  possible  eigenvalues  were  statistically  significant.  The  per- 
turbation reduces  the  overall  rms  error  by  only  3 0 percent;  thus,  it  is  obvious  that  there  also 
must  exist  velocity  variations  of  a shorter  wavelength.  The  rms  variation  of  perturbation  in  each 
shell  is  small  — it  ranges  from  0.3  percent  for  the  upper  mantle  to  0.1  percent  for  Region  III.  The 
fact  that  lateral  variations  are  the  least  in  the  middle  mantle  is  in  agreement  with  the  results  of 

4 

Julian  and  Sengupta,  who  used  only  deep  events;  the  data  in  this  study  were  dominated  by  shallow 
focus  sources. 

Figures  IV-1  (a)  through  (d)  show  lateral  variations  in  velocity  obtained  by  the  least-squares 
fit  of  spherical  harmonics  up  to  the  third  degree  to  results  in  Table  IV-1;  the  areas  of  individual 
blocks  were  used  as  weighting  factors.  It  is  interesting  to  notice  that  velocity  variations  in 
Regions  III  and  IV  are  strongly  correlated.  The  zones  in  the  east  Pacific  and  southern  Asia 
are  anomalous  at  all  levels. 

Several  experiments  have  been  performed  to  check  the  stability  of  the  solution.  First, 
results  were  compared  using  data  with  and  without  corrections  for  ellipticity.  In  the  latter  case, 
one  obtains  velocity  perturbations  strongly  influenced  by  the  zonal  harmonic  of  degree  2 — exactly 
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TABLE  IV-1 

RESULTS  OF  INVERSION  FOR  LATERAL  VARIATIONS  IN  COMPRESSIONAL  VELOCITY.  AVERAGE  PERTURBATION  OF  VELOCITY  AV 
IS  REMOVED  FROM  THE  VALUES  REPORTED  FOR  INDIVIDUAL  BLOCKS. 
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as  expected.  To  check  whether  the  large  differences  in  the  number  of  rays  sampling  individual 
blocks  may  introduce  a bias,  the  data  were  counterweighted  according  to  the  number  of  connec- 
tions between  the  blocks  containing  the  source  and  the  receiver.  The  resulting  change  in  the 
pattern  of  velocity  anomalies  was  insignificant.  An  allowance  also  was  made  for  the  possibility 
that  the  radius  of  the  670-km  discontinuity  may  be  variable.  In  this  case,  the  change  in  radius 
represented  an  additional  30  parameters.  The  number  of  significant  eigenvalues  was  132  out  of 
150  possible.  This  indicates  that  the  resolution  for  this  parameter  is  relatively  poor,  however, 
only  results  for  Region  II,  immediately  underlaying  the  discontinuity,  were  appreciably  changed. 
Perturbations  in  the  radius  occasionally  reach  40  km;  these  are  rather  large,  but  differences 

5 

of  tens  of  kilometers  have  been  reported  on  the  basis  of  observations  of  P'P'  reflections. 

Finally,  Hager^  took  the  results  for  the  case  with  a variable  radius  of  the  670-km  discon- 
tinuity, assumed  that  density  perturbations  are  proportional  to  velocity  variations  according  to 

7 8 

Birch's  law,  computed  equivalent  geoid  heights,  and  compared  them  with  the  observed  heights. 

He  has  found  that  there  are  statistically  significant  (at  the  98-percent  level)  correlations  between 
the  observed  geoid  heights  (with  respect  to  the  hydrostatic  spheroid  of  reference)  and  separate 
contributions  from  the  upper  and  lower  mantles.  These  contributions  are  of  opposite  sign,  such 
that  approximately  80  percent  of  the  effect  is  canceled  by  addition.  The  sum  of  the  contributions 
does  not  show  significant  correlation  with  observations  because  of  the  noise  in  the  velocity  struc- 
ture, but  it  has  the  correct  magnitude. 

These  results,  if  confirmed  by  further  tests,  could  be  of  fundamental  importance.  Gravita- 
tional compensation  between  the  upper  and  lower  mantles  indicates  transport  of  the  material 
between  these  two  regions.  Good  correlation  with  the  gravity  anomalies,  calculated  with  respect 
to  the  hydrostatic  ellipsoid  rather  than  the  ellipsoid  of  observed  flattening,  indicates  that  the 
"equatorial  bulge"  is  an  expression  of  a lateral  heterogeneity  rather  than  a residuum  of  the 

earth's  flattening  when  it  rotated  with  higher  angular  velocity.  The  latter  interpretation  has  led 
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to  estimates  of  the  lower-mantle  viscosity  as  high  as  6x10  stokes  and  was  used  as  a strong 

9 

argument  against  convection  in  the  lower  mantle. 

Even  though  the  results  presented  here  are  encouraging,  it  is  necessary  to  perform  further 
experiments  before  firm  conclusions  may  be  drawn.  In  particular,  tests  are  needed  to  investi- 
gate the  stability  of  the  solution  with  respect  to  major  changes  in  the  discretization  system. 

A.  M.  Dziewonski 

B.  DETERMINATION  OF  THE  THREE-DIMENSIONAL  STRUCTURE 
OF  THE  LITHOSPHERE 

A new  approach  to  three-dimensional  seismic  modeling  of  the  lithosphere  has  been  studied. 

In  contrast  to  conventional  methods,  the  assumption  of  uniform  material  property  within  a layer 
is  relaxed,  making  it  possible  to  accommodate  more  complex  geological  structures.  Another 
important  feature  is  that  we  can  estimate  a lower  limit  of  the  true  rms  velocity  perturbations 
in  the  lithosphere,  under  the  assumption  of  ray-theory.  Using  NORSAR  P-wave  residuals  for 
teleseismic  events,  we  have  estimated  the  three-dimensional  seismic  anomalies  in  the  litho- 
sphere to  the  depth  of  126  km  beneath  the  array.  The  fit  to  the  observations  is  excellent,  as 
indicated  in  the  final  residual  comparable  to  the  measurement  error.  We  found  clear  evidence 
of  pipe -like  structures  below  Moho,  dipping  northward  and  away  from  the  surface  contours  of 
the  Oslo  Graben  (see  Fig.  IV-2).  These  pipe-like  structures  were  interpreted  as  vestiges  of 
magma  ascent  due  to  penetrative  convections  associated  with  the  Permian  volcanism  of 
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Oslo  Graben.  The  northward  inclination  of  the  pipe-like  structures  is  interpreted  as  a result 
of  plastic  deformation  of  the  lithosphere  due  to  the  shear  exerted  by  the  asthenosphere  convection 
current  driving  the  plate  motion. 

Our  results  also  show  that  the  velocity  fluctuations  are  of  the  same  order  in  the  crust  as  in 
the  upper  mantle,  at  least  to  the  depth  of  126  km.  The  rms  of  the  true  velocity  fluctuations  in 
the  crust  and  upper  mantle  under  NORSAR  is  at  least  3.4  percent,  and  is  in  agreement  with 
estimates  obtained  from  a statistical  analysis  of  P-time  fluctuation  based  on  the  Chernov  theory. 

K.  Aki 

L.  A.  Christoffersson 

E.  Husebye 


C.  THREE-COMPONENT  ANALYSIS  OF  SEISMIC  CODAS 
FROM  NOVAYA  ZEMLYA  EXPLOSIONS 

Six  presumed  nuclear  explosions  from  Novaya  Zemlya  were  studied  using  three-component, 
short-period  data  from  the  D2  subarray  at  LASA.  The  results  show  that  each  coda  contains  a 
large  number  of  strong,  impulsive  body  phases  scattered  by  inhomogeneities  in  and  out  of  the 
diametral  plane.  In  particular,  strong  phases  on  the  SH  and  SV  components  arrive  several 
seconds  after  the  P-arrival,  and  continue  to  arrive  for  at  least  4 minutes  of  coda  for  each  event. 

Two  northern  events,  with  body  wave  magnitudes  6.7  and  6.8,  have  remarkably  similar 
codas,  suggesting  identical  site  locations  for  the  two  events.  The  remaining  two  northern 
events,  with  magnitudes  6.1  and  6.3,  have  codas  which  match  each  other  reasonably  well,  but 
do  not  match  the  codas  of  the  first  two  events.  To  the  south,  the  two  events  have  magnitudes 
6.0  and  6.9,  the  codas  which  do  not  clearly  match  each  other  or  any  of  the  codas  of  the  northern 
events.  This  indicates  that  small  changes  in  source  location  or  site  characteristics  may  com- 
pletely alter  the  complexity  of  the  coda  arriving  at  distant  stations. 

Table  IV-2  is  a list  of  the  six  events,  and  Fig.IV-3  shows  their  PDE  locations  as  black  dots. 
Four  of  the  events  are  clustered  at  approximately  73  °N  55 °E,  and  the  remaining  two  are  located 
further  south  near  71  °N  54 °E.  The  three -component  data  were  rotated  to  radial  horizontal  (R), 

vertical  (Z),  and  tangential  horizontal  (T)  directions,  and  filtered  to  enhance  polarized  body 
10 

waves. 


TABLE  IV-2 

LIST  OF  PRESUMED  NUCLEAR  EXPLOSIONS 
AT  NOVAYA  ZEMLYA 

Origin  Time 

Latitude 

Longitude 

Event 

Da  te 

(GMT) 

(°N) 

(°E) 

1 

14  Oct  1969 

7:00:06 

73.40 

54.81 

2 

28  Aug  1972 

5:59:57 

73.34 

55.09 

3 

14  Oct  1970 

5:59:57 

73.32 

55. 15 

4 

12  Sep  1973 

6:57:43 

73.30 

55.16 

5 

27  Sep  1973 

6:59:58 

70.76 

53.87 

6 

27  Oct  1973 

6:59:57 

70. 78 

54.18 
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Figure  IV-4  shows  the  three  components  of  Event  1 before  and  after  the  polarization  filtering. 
A striking  feature  of  the  unfiltered  data  is  the  high  levels  of  coda  energy  on  the  R and  T compo- 
nents minutes  after  the  P-wave.  The  filtered  components  show  a succession  of  impulsive  events 
arriving  throughout  the  coda.  These  events  are  highly  polarized  P-  and  S-body  phases  of  all 
azimuths  and  dips.  The  suppressed  portions  of  the  coda  cannot  be  dominantly  P-energy  radiated 
from  the  source  region  but  must  be  overlapping  arrivals  of  scattered  P-  and  S-waves  from  dif- 
ferent regions  as  well  as  locally  generated  surface  waves,  which  produce  particle  motions  which 

are  not  clearly  polarized.  This  suggests  that  scattering  regions  are  more  numerous  and  complex 

11 

than  usually  deduced  from  array  studies  of  Z-components  only. 

Figures  IV- 5,  IV-6,  and  IV-7  display  the  polarized  components  Z,  R,  and  T of  all  six 
events,  each  scaled  to  m^  = 6.5.  In  Fig.  IV- 5,  the  Z-components  of  Events  3 and  4 with  = 6.7 
and  6.8,  respectively,  are  nearly  identical  for  the  duration  of  the  coda  displayed.  This  suggests 
that  the  two  events  have  close  or  identical  site  locations,  as  indicated  by  the  PDE  locations  in 
Fig.  IV-3.  The  Z-components  of  Events  1 and  2,  with  magnitudes  6.1  and  6.3,  respectively,  also 
match  but  not  as  well  as  Events  3 and  4. 

The  four  events  were  relocated  using  a common  set  of  nineteen  WWSSN  stations.  The  re- 
vised locations  are  shown  by  circles  in  Fig.  IV-3.  It  appears  that  the  events  group  themselves 
into  two  pairs  which  are  separated  by  10  km.  Although  the  absolute  locations  may  be  question- 
able, the  separation  between  pairs  is  probably  real.  The  difference  in  arrival  times  between 
the  pairs  is  on  the  average  about  3-  to  10-times  larger  than  the  variation  in  station  residuals 
at  the  stations  most  sensitive  to  the  geographic  separation. 

The  mismatch  in  polarized  codas  between  the  two  pairs  of  events  suggests  that  10  km  is  a 
large  enough  separation  to  cause  the  codas  to  be  uncorrelated.  Correlation  coefficients  between 
similar  components  of  Events  3 and  4 are  very  high,  typically  >0.8  for  4 minutes  of  data  using 
windows  50  sec  long.  Events  1 and  3,  however,  have  essentially  uncorrelated  components  with 
correlation  coefficients  between  —0.3  to  +0.3. 

Events  5 and  6,  located  at  the  southern  end  of  Novaya  Zemlya,  are  separated  by  10  km 
according  to  the  PDE  locations.  Using  a set  of  20  stations,  they  were  both  relocated  about  10  km 
to  the  west  with  no  change  in  separation.  Figure  IV-5  shows  that  the  polarized  codas  of  these 
two  events  do  not  match  each  other  or  the  codas  of  the  northern  events.  This  suggests  that 
Events  5 and  6 are  separated  by  more  than  the  correlation  distance  for  scattered  waves  produced 
by  the  heterogeneities  near  Novaya  Zemlya. 

Beams  of  Events  1 to  4 are  displayed  in  Figs.  IV-5,  IV-6,  and  IV-7.  Some  of  the  phases 
shown  on  the  beams  are  very  coherent  for  all  four  events.  These  phases  represent  strong 
scatterers  in  the  mantle,  and  are  being  studied  in  detail. 

C.  W.  Frasier 
M.  Yang 
R.  E.  Needham 

D.  SCATTERING  OF  ELASTIC  WAVES  FROM  DEPTH-DEPENDENT 
INHOMOGENEITIES 

A new  approximate  solution  has  been  developed  for  calculating  the  reflection  and  transmission 
coefficients  for  plane  elastic  waves  incident  to  a transition  zone  in  which  elastic  parameters  vary 
with  depth.  The  method  is  approximate  in  that  primary  scattering  with  no  multiples  is  calculated 
at  each  depth  point  of  the  transition  zone.  The  effect  of  the  total  transition  zone  is  obtained  by 
integrating  the  primary  scattered  field  over  the  total  depth  thickness  of  the  zone.  Surprisingly, 


39 


the  method  seems  to  be  accurate  for  transition  zones  in  which  elastic  parameters  vary  smoothly 
except  at  frequencies  near  0 Hz. 

Briefly,  the  method  associates  with  each  element  AZ  in  the  transition  zone  a scattered  wave 

which  may  contribute  to  either  a reflection  or  transmission  response.  For  the  reflection  re- 

1 2 

sponse  rpp,  for  example,  the  scattered  amplitude  from  an  element  AZ  is  given  by  Bortfeld  as 


Arpp  = 


cos  (2j) 


ii  Bl 

2>  p 


; — z- 

2 cos  l 


SL  - 4 
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sin  j 


P 1 

T 
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(IV-1) 


where  i and  j are  incident  P-  and  SV-wave  angles,  respectively,  at  depth  Z,  and  primes  (') 
indicate  derivatives  of  parameters  with  respect  to  Z.  The  phase  factor  associated  with  this 
scattered  amplitude  is  given  by 


exp[  — ico(sx  + 2t  ^)] 


(IV-2) 


where  s is  the  horizontal  slowness  of  the  wave  in  the  x-direction,  and  t is  the  one-way  vertical 
delay  time  from  the  top  of  the  transition  zone  to  the  scattering  element  AZ.  The  total  reflected 
primary  response  is  obtained  by  integrating  the  product  of  Eqs.  (IV-1)  and  (IV-2)  over  the  tran- 
sition zone  thickness.  Thus, 


RPP(co) 


— (rpp)  exp[  — icj(sx  + 2^)]  dZ 


(IV-3) 


By  a change  of  integration  variable  from  depth  Z to  arrival  time  t = sx  + 2t(  , Eq.  (IV-3)  can  be 
written  as 


RPPM  = £ (rpp)  e'iwt  dt  . (IV-4) 

^ - oo 

This  equation  can  be  inverted  directly  to  time  by  inspection,  yielding 

PPP(t)  = -pL  (rpp)  ^ — r 
dZ  2 cos  l 

where  Z can  be  calculated  as  a parametric  function  of  t by  evaluating  the  expression 

1 - “ ■ 2 1,Z  dz  • t'v-6> 
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Similar  formulas  can  be  found  for  the  other  primary  scattered  responses  RPS(t),  RSS(t), 
RSP(t),  TPS(t),  and  TSP(t).  It  can  be  shown  that  the  responses  TPP(t)  and  TSS(t)  contain  no 
primary  interactions,  so  that  they  are  approximated  by  weighted  delta  functions. 

To  illustrate  this  method,  we  computed  reflection  and  transmission  responses  for  three 
transition- zone  models  which  are  tabulated  in  Table  IV-3.  In  Model  1,  the  density  varies 
linearly  with  depth;  Model  2 contains  shear  and  compressional  velocity  gradients;  and  Model  3 
has  gradients  in  density  and  each  velocity.  In  each  model,  the  transition  zone  is  5 km  thick 
between  two  elastic  half-spaces. 

Figures  IV -8,  IV-9,  and  IV-10  show  the  primary  pulse  shapes  for  Models  1,  2,  and  3, 

respectively.  Also  displayed  in  these  figures  are  the  total  reflection  and  transmission  response 

1 3 

shapes  calculated  by  Haskell's  layer  matrix  method.  From  Fig.  IV-8,  it  is  clear  that  the 


z=zm 


( I V - 5 ) 
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PARAMETERS 

TABLE  IV-3 

OF  UPPER  AND  LOWER  HALF-SPACES 
FOR  EACH  MODEL 

Model  1 

Model  2 

Model  3 

o 

Q_ 

2.0 

3.0 

2.0 

Upper  half-space  - 

ao 

5.0 

5.0 

5.0 

A J 

3.0 

3.0 

3.0 

Transition  zone  5 km  thick 

'pi 

3.5 

3.0 

3.5 

Lower  half-space  4 

ai 

5.0 

8.0 

8.0 

.Pi 

3.0 

4.8 

4.8 

primary  scattered  field  for  Model  1,  containing  only  a density  gradient,  is  a very  good  approxi- 
mation to  the  total  scattered  field.  Figures  IV -9  and  IV-10  show  some  discrepancy  between 
primary  and  total  pulse  shapes,  the  largest  mismatches  being  for  Model  3 in  Fig.  IV-10.  The 
critical  incident  P-wave  angle  for  Models  2 and  3 is  38.7°. 

In  order  to  investigate  the  higher-order  scattering  neglected  in  the  primary  pulse  calcula- 

14 

tions,  we  computed  the  responses  for  Model  3 using  a discrete  time  method.  This  method, 
like  Haskell's  technique,  requires  a set  of  thin  homogeneous  layers  to  approximate  a transition 
zone.  By  adjusting  the  layer  thicknesses,  the  primary  arrivals  were  constrained  to  arrive  at 
equal  time  intervals  as  shown  by  the  tall  spikes  in  Figs.IV-11,  IV-12,  and  IV-13.  Between  the 
large  spikes  are  higher-order  scattered  arrivals  which,  although  low  in  amplitude,  can  be  dense 
in  time.  The  result  of  this  is  that  the  total  field,  obtained  by  convolving  the  discrete  time  re- 
sponses with  a smoothing  kernel,  may  look  different  from  the  primary  pulse  shape.  This  is 
most  evident  for  the  reflection  responses  at  non-normal  incidence,  e.g.,  RSP  and  RSS. 

Figures  IV-14  and  IV -15  show  the  amplitude  responses  for  Model  3 from  0 to  2.5  Hz,  calcu- 
lated by  three  methods  — two  of  which  are  the  Haskell  layer  matrix  method  and  the  primary  pulse 
method  of  this  report.  These  responses  are  Fourier  transforms  of  the  pulse  shapes  displayed 
in  Fig.  IV-10. 

Also  given  is  the  0(l/u>)  response,  equal  to  the  spectrum  of  the  initial  and  final  discontinuities 
of  the  primary  pulse  shapes  shown  in  Fig.  IV-10.  Except  for  RSS  at  non-normal  incidence,  the 
spectrum  of  the  discontinuities  comprises  the  major  component  of  the  total  scattered  spectrum. 

For  most  short-period  instruments,  the  ground  velocity  response  rapidly  falls  off  at  0.5  Hz 
and  below.  This  implies  that  the  primary  pulse  method  described  in  this  report  is  a quick  and 
accurate  way  to  compute  reflection  and  transmission  responses  for  sharp  transition  zones  in  the 
earth  recorded  on  conventional  instruments. 

C.  W.  Frasier 
P.  Richards^ 


t Lamont-Doherty  Geological  Observatory,  Palisades,  New  York  10964. 
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E.  GLOBAL  SCATTERING  PARAMETERS  IN  THE  SHORT -PERIOD  BAND 


15  10 

Previous  work  f has  demonstrated  that  the  short-period  codas  of  both  local  and  tele- 
seismic  events  observed  at  LASA  contain  considerable  energy  attributable  to  scattering  by 
localized  inhomogeneities.  In  order  to  formulate  an  improved  model  for  scattering,  it  is  nec- 
essary to  determine  the  limits  of  validity  of  currently  used  models;  specifically,  the  application 
of  Chernov’s  work  on  acoustic  waves  in  random  liquid  media  to  the  study  of  seismic  scattering 
must  be  examined.  To  do  this,  a program  has  been  undertaken  to  determine  average  values  for 

the  earth  of  the  scattering  parameters  used  in  Chernov's  work. 

1 6 

Aki  has  shown  that  for  waves  with  frequencies  greater  than  1 Hz,  Chernov’s  work  breaks 
down  when  applied  to  LASA;  so,  in  this  study,  work  is  being  restricted  to  P-waves  with  periods 
of  0.25  to  1 sec.  Data  from  LASA,  NORSAR,  and  Russian  stations  will  be  used.  The  configura- 
tion of  LASA  and  NORSAR  allows  the  study  of  various  scales  of  inhomogeneities.  The  very  fine 
structure  is  examined  by  using  seismometers  in  the  A and  B rings  with  a mean  spacing  of 
15  km.  Coarser  structure  is  examined  by  using  subarray  sums  from  the  A,  C,  D,  E,  and  F 
rings  with  an  effective  mean  spacing  of  50  km. 

The  first  15  events  studied  cover  the  region  along  a rough  arc  from  Alaska  through  the 
Aleutians,  Kamchatka,  Japan,  the  Philippines,  Borneo,  Sumatra,  and  the  Banda  Sea  to  the  west 
to  the  Red  Sea.  These  events  will  be  used  for  the  initial  determination  of  the  mean  residuals 
near  the  receivers,  and  can  be  used  for  the  study  of  shallow  inhomogeneities  near  the  receivers. 

The  next  15  events  were  selected  from  the  Kurile  Islands  and  Japan  region.  One -half  of 
these  events  are  shallow  (<60  km  deep),  and  the  other  half  are  deep  (>150  km).  The  events  were 
chosen  along  a great  circle  to  LASA  and  perpendicular  to  the  great  circle  to  NORSAR.  All  the 
events  are  from  65°  to  90°  from  LASA,  to  avoid  the  core  shadow,  and  about  80°  from  NORSAR. 
These  events  are  at  least  magnitude  5.5  and  exhibit  a clean  impulsive  P-arrival.  They  will  be 
processed  using  methods  developed  in  concert  with  Dr.  A.  V.  Nikolaev  of  the  Institute  of  Physics 
of  the  Earth  and  will  be  used  to  study  variations  in  the  amplitude,  travel-time,  and  waveform 
of  the  P-arrival.  These  events  will  yield  information  on  inhomogeneities  which  lie  in  the  source 
region. 

The  final  step  in  the  study  will  be  to  attempt  to  study  lower-mantle  inhomogeneities.  For 
this  purpose,  15  events  will  be  selected  along  a line  from  Japan  to  the  southwest  toward  Taiwan 
and  the  Philippines.  Rays  from  earthquakes  with  magnitude  >6.5  diffracted  in  the  D’  or  D"  layers 
should  give  information  as  to  whether  the  degree  of  inhomogeneity  changes  between  these  two 
layers.  The  rate  at  which  the  amplitude  of  PkP  falls  off  after  90°  can  be  used  to  determine 
whether  the  region  near  to  the  core-mantle  boundary  is  more  strongly  inhomogeneous  than  the 
rest  of  the  mantle. 

J.  Scheimer 
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(a)  Region  I 


(b)  Region  II 


Fig.  IV-1  (a-d).  Lateral  variations  of  compressional  velocities  in  four  concentric 
shells;  contour  units  are  in  meters/second.  Depth  range  of  individual  regions  is: 
Region  I,  from  0 to  670  km;  Region  n,  from  670  to  1400  km;  Region  III,  from 
1400  to  2200  km;  and  Region  IV,  from  2200  km  to  core-mantle  boundary. 
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(c)  Region  III 


(d)  Region  IV 


Fig.  IV-1  (a-d).  Continued. 
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10  km  SOUTH  OF  AO 


Fig.  IV-2.  NORSAR  SW-NE  cross  section  10  km  south  of  AO. 


3 14  OCT  1970  6.7 

4 12  SEP  1973  6.8 

5 27  SEP  1973  6.0 

6 27  OCT  1973  6.9 


Fig.IV-3.  Map  of  Novaya  Zemlya  showing  locations  of  six  presumed 
explosions.  Black  circles  are  PDE  locations,  and  open  circles  are  four 
relocated  epicenters  using  19  WWSSN  stations. 
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Fig.  IV-4.  Three-component  data  for  Event  1.  Above  unfiltered  data 
is  polarization  gain  which  equals  1 if  data  are  linearly  polarized,  and 
0 if  circularly  polarized.  Polarized  data  equal  principal  component  of 
data  projected  onto  data  and  scaled  by  polarization  gain. 
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Fig.IV-5.  Z components  of  six  events,  each  polarized  and  normalized 
to  = 6.5.  Also  shown  is  sum  of  Events  1 to  4,  sum  of  Events  5 
and  6,  and  sum  of  Events  1 to  6. 
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IV-6.  R components  of  six  events,  each  polarized  and  normalized  to  m^  = 6.5. 
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IV-7.  T components  of  six  events,  each  polarized  and  normalized  to  = 6.5. 
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MODEL  1 


HASKELL  PRIMARY  PULSE 
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Fig.  IV-8.  Time  domain  reflection  and  transmission  responses  for  Model  1.  Each 
response  is  computed  by  Haskell1  s method  and  primary-pulse  method,  which  are 
displayed  side  by  side.  All  angles  are  incident  P-wave  angles.  SV  incident  waves 
have  same  slowness  as  given  P-angle  in  each  response. 
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MODEL  2 


HASKELL  PRIMARY  PULSE  HASKELL  PRIMARY  PULSE 


Fig.  IV-9.  Time  domain  reflection  and  transmission  responses  for  IVIodel  2. 
As  in  Fig.  IV -8,  all  angles  displayed  are  P-wave  incident  angles. 
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Fig.  IV-10.  Time  domain  reflection  and  transmission  responses  for  Model  3.  All 
angles  are  incident  P-wave  angles. 


51 


PRIMARY  DISCRETE 

PULSE  TIME 


DISCRETE  HASKELL 

TIME 

(smoothed) 


L 

o 


(sec) 


iw-Hw-n 


Fig.  IV-11.  Time  domain  reflection  responses  RPP  for  Model  3 calculated 
by  three  different  methods  for  normal  incidence  and  near  critical  incidence. 
Smoothed  discrete  time  method  is  obtained  by  convolving  unsmoothed  re- 
sponse with  a triangular  kernel  of  area  1. 
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Fig.  IV-12.  Time  domain  reflection  responses  RPS  and  RSP  for  Model  3 
at  near  critical  incident  angles.  Angles  given  are  P-wave. 
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Fig.  IV-13.  Time  domain  reflection  response  RSS  for  Model  3 at  normal 
and  near  critical  incident  angles.  Angles  given  are  for  P-waves  of  same 
slowness  as  incident  SV-waves. 
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Fig.  IV-14.  Frequency  domain  responses 
RPP  and  RSS  for  Model  3 at  normal  and 
near  critical  angles  of  incidence.  Each 
0 ( l/c o)  response  is  spectrum  of  time  do- 
main discontinuities  at  beginning  and  end 
of  primary  pulse  response. 


Hz 

Fig.  IV-15.  Frequency  domain  responses 
RPS,  RSP,  TPS,  and  TSP  at  near  critical 
angle  of  incidence  to  Model  3.  Each  0(l/o;) 
response  is  spectrum  of  two  time  domain 
discontinuities  in  primary  pulse  shape. 
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V.  EARTHQUAKE  OCCURRENCE  RATES 


A.  SPATIAL  AND  TEMPORAL  VARIATIONS 

IN  THE  FREQUENCY-MAGNITUDE  CURVE 

Spatial  variations  in  the  level  of  the  frequency-magnitude  curve  are  well  known,  and  in 
Secs.  B and  C below,  the  temporal  variations  in  this  level  are  discussed.  First,  however,  we 
shall  examine  evidence  for  changes  in  the  shape  of  the  curve  in  space  and  time.  In  spite  of  a 
voluminous  literature,  the  evidence  that  these  changes  in  shape  occur  is  still  very  questionable. 

The  detection  of  spatial  variations  in  shape  is  confused  by  differences  in  the  definition  of 
magnitude,  and  by  network  detection  characteristics  (which  are  a strong  function  of  location  at 
lower  magnitudes).  If  they  are  real,  spatial  variations  should  be  readily  apparent  in  the  NEIS 
PDE  catalog,  which  has  a good  global  detection  threshold  and  a reasonably  consistent  measure 
of  magnitude  m^. 

Data  from  this  catalog  for  the  period  1964.0  to  1974.3  have  been  subdivided  on  the  basis  of 
seismic  region  number  into  three  broad  tectonic  types:  subduction  zones,  oceanic  ridge  and 
transform  areas,  and  shield  areas.  Frequency-magnitude  curves  for  these  three  regimes  arc 
shown  in  Fig.  V-l.  In  spite  of  the  large  differences  in  level,  these  curves  match  one  another 
very  closely,  except  for  inevitable  scatter  at  high  magnitudes  where  there  are  very  few  events. 

A representation  of  these  types  of  data  that  recently  has  been  introduced  is  the  b-m.  plot. 

b ^ 

The  maximum-likelihood  estimate  of  the  slope  b of  a straight  line  fitted  to  the  data  is  given  by 
0.4343 


b = 


M - M 


where  M is  the  average  magnitude  of  all  the  events  considered,  and  M^.^  is  the  lower  thresh- 
old of  magnitude  used.  Clearly,  an  estimate  of  b can  be  made  for  various  different  thresholds, 
and  then  b can  be  plotted  against  threshold. 

Figure  V-2  shows  such  plots  for  the  three  tectonic  environments.  95-percent  confidence 
limits  on  each  point  may  be  calculated  as  ±1.96  b/\rn,  where  n is  the  number  of  events  used. 
The  differences  between  these  curves  are  not  significant  at  the  95-percent  level.  Notice,  how- 
ever, that  all  three  depart  significantly  from  linearity. 

We  have  concluded  from  these  results,  and  from  similar  data  from  smaller  regions,  that 
the  10  years  of  data  in  the  PDE  listing  do  not  indicate  significant  variations  in  space  of  the  Shape 
of  the  frequency-magnitude  curve.  Studies  in  the  literature  that  appear  to  show  otherwise  (e.g. , 
Acharya  ) may  be  entirely  the  result  of  failing  to  take  account  of  the  spatial  variation  in  net- 
work detection  threshold. 

Temporal  variations  in  the  shape  of  the  frequency-magnitude  curve  have  received  less 
attention  in  the  literature.  We  show  one  example  from  the  PDE  catalog  that  argues  against  these 
variations  on  a time  scale  of  years.  Figure  V-3  shows  the  frequency-magnitude  curve  for  the 
Northwest  Pacific  during  the  activity  maximum  of  1964-65,  and  for  the  period  1966-74,  both 
normalized  to  events  per  year.  The  differences  between  the  shape  of  these  curves  arc  not  sig- 
nificant. This  result  makes  an  interesting  comparison  with  those  described  in  Sec.  B below. 

M.  A.  Chinnery 


55 


B.  CORRELATIONS  BETWEEN  SEISMIC  ACTIVITY 
IN  WIDELY  SEPARATED  REGIONS 


,4  . 


The  Lincoln  Laboratory  Data  Analysis  and  Display  System  (DADS)  is  a powerful  tool  for  the 
examination  of  tabular  data.  With  the  aid  of  this  system,  we  have  carried  out  a number  of  inves- 
tigations of  the  characteristics  of  earthquake  catalogs.  The  present  study  compares  temporal 
variations  in  the  level  of  seismic  activity  between  different  areas. 

The  principal  source  of  data  for  this  study  is  the  NEIS  PDE  catalog.  This  listing  appears 
to  have  been  relatively  homogeneous  since  1964,  and  we  have  therefore  used  this  list  for  the 
period  1 January  1964  to  30  April  1974.  Supplementing  this,  we  have  the  Regional  Catalogue  of 
Earthquakes  published  by  the  Japan  Meteorological  Agency  (January  1964  to  January  1970),  the 
list  of  earthquakes  in  the  Southern  California  region  published  by  the  California  Institute  of 
Technology  (January  1964  to  December  1971),  and  the  Bulletin  of  the  International  Seismological 
Center  (January  1964  to  December  1970). 

These  catalogs  were  used  to  construct  earthquake  time  series  for  a number  of  regions,  se- 
lected on  the  basis  of  seismic  region  number.  These  time  series  consist  of  counts  of  earthquakes 
per  unit  time  with  magnitudes  greater  than  or  equal  to  a given  threshold.  Thresholds  were  chosen 
so  that  reasonably  complete  counts  could  be  expected.  In  the  case  of  the  PDE  listing,  earthquakes 
smaller  than  m^  = 5.2  were  not  used. 

Data  shown  in  the  present  study  are  obtained  from  counts  over  40 -day  intervals.  An  example 
of  one  such  histogram  is  shown  in  Fig.  V-4(a)  for  shallow  subduction  zone  earthquakes  bordering 
the  Pacific  Ocean.  The  time  series  constructed  in  this  way  contain  a considerable  amount  of  high- 
frequency  "noise,”  and  the  character  of  the  time  series  becomes  clearer  after  smoothing.  We 
have  applied  a sliding  200-day  window  to  the  data,  and  the  result  is  shown  in  Fig.  V-4(b). 

The  activity  maximum  in  1965  is  only  partially  the  result  of  the  Alaskan  earthquake  and  the 
Rat  Island  sequence.  A similar  graph  [Fig.  V-4(c)]  for  deep  earthquakes  (depth  > 300  km)  also 
shows  a maximum  during  1965.  In  both  cases,  the  maximum  is  a factor  of  4 or  5 higher  than 
the  subsequent  lows  in  mid-1966  and  mid-1967,  and  a factor  of  about  2 higher  than  the  mean  for 
the  whole  period. 

, Time  series  have  been  constructed  for  a number  of  smaller  regions  of  the  globe,  and  12  of 
these  are  shown  in  Figs.  V-5  and  V-6.  All  regions  show  an  activity  maximum  duringl965,  followed 
by  minima  in  mid-1966  and  mid-1967.  The  location  of  the  1965  peak  varies  considerably  with 
region,  but  this  may  be  a local  phenomenon,  since  the  subsequent  minima  show  much  less  vari- 
ation. After  1968,  the  correlation  between  the  time  series  is  much  less  impressive. 

These  data  are  interpreted  as  a first  indication  of  the  presence  of  a global  triggering  stress 
that  fluctuates  rapidly.  Stress  release  by  large  earthquakes  will  normally  occur  near  peaks  of 
the  triggering  stress,  when  a region  is  very  close  to  failure.  When  the  local  stress  is  not  close 
to  failure,  variations  in  the  applied  stress  appear  to  be  reflected  in  the  numbers  of  small-  to 
medium-size  earthquakes  that  occur.  This  suggests  that  information  about  the  nature  of  the  trig- 
gering stress  may  be  obtained  from  earthquake  time  series. 

Detailed  study  of  these  time  series  is  continuing.  Their  spectral  characteristics  are  dis- 

cussed  in  Sec.  C below.  M.  A.  chinnery 

T.  E.  Landers 
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C.  SPECTRA  FROM  TWO  MONTHS  TO  TWO  YEARS  OF  REGIONAL 
EARTHQUAKE  OCCURRENCE  RATES 

In  a previous  SATS^  the  level  of  seismic  activity  (events  per  magnitude  interval)  was  shown 
to  be  a time -dependent  function  independent  of  the  magnitude  range  over  which  it  is  determined. 
Consequently,  the  spectrum  of  the  corresponding  seismic  occurrence  rate  function  (events  within 
a magnitude  range  per  time  interval)  for  global  seismicity  data  was  examined  and  found  to  have 
a significantly  nonrandom  concentration  of  power  at  a period  near  230  days.7  The  energy  con- 
tained in  this  spectral  peak  amounted  to  approximately  one -third  of  the  total  energy  in  the  spec- 
trum. To  determine  whether  this  periodicity  is  a global  phenomenon  of  all  seismic  regions  or 
peculiar  to  certain  tectonic  regions,  the  spectra  of  the  earthquake  occurrence  rate  sequences  for 

various  tectonic  areas  were  computed.  In  addition,  the  resolution  of  the  maximum  entropy  spcc- 
7 

tral  analysis  was  increased  to  determine  if  other  less-powerful  nonrandom  spectral  peaks  were 
present. 

The  following  is  a list  of  zones  considered. 

o 

(1)  The  entire  shallow  (depth  100  km)  Pacific  (Regions  1 to  301). 

(2)  The  entire  nonshallow  Pacific. 

(3)  The  entire  shallow  Atlantic  (Regions  402  to  414). 

(4)  Shallow  India-Tibet-Szechwan- Yunan  (Regions  302  to  319). 

(5)  Shallow  Western  Asia-Middle  East-Crimea-Balkans-Western  Medi- 
terranean (Regions  347  to  401). 

(6)  The  following  Pacific  zones,  all  shallow: 

(a)  Alaska-Aleutian  Arc  (Regions  1 to  17) 

(b)  Mexico-Guatemala  Area-Central  America-Caribbean 
Loop-Andean  South  America  (Regions  53  to  146) 

(c)  Kcrmadec-Tonga-Samoa  Area -Fiji  Islands  Area 
(Regions  169  to  182) 

(d)  New  Hebrides  Islands-Bismark  and  Solomon  Islands- 
New  Guinea-Caroline  Islands  to  Guam  (Regions  183 
to  210). 

(e)  Japan-Kuriles -Kamchatka  (Regions  217  to  230) 

(f)  Philippines -Borneo-Celebes-Sunda  Arc-Burma- 
Southeast  Asia  (Regions  248  to  301). 

The  data  source  was  the  PDE,  and  the  time  span  was  1964  to  the  present. 

The  sequences  were  prepared  in  the  following  manner.  To  prevent  aliasing,  daily  sum  data 
were  lowpassed  with  a 3-pole  Butterworth  filter  with  a 3-dB  cutoff  at  l/(40  days),  and  then  sam- 
pled at  20-  or  10-day  intervals.  To  preserve  the  phase,  the  filter  was  passed  first  forward  and 
then  reversed  and  the  data  truncated  at  both  ends  at  twice  the  impulse  response  of  the  filter.  In 
a similar  manner,  the  data  were  highpassed  above  l/(2  years).  The  resulting  sequences  were 

9 

then  analyzed  using  the  Burg  Maximum  Entropy  Method  (BMEM). 
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The  most-critical  and  least -understood  parameter  in  computing  a spectrum  using  the  BMEM 
is  the  filter  length.  The  use  of  a filter  with  a length  greater  than  the  order  of  the  process  can 
result  in  spurious  peaks.  A shorter  filter  can  diminish  the  resolution.  Filter  lengths  comparable 

Q 

to  the  number  of  data  points  for  high-order  processes  can  bias  the  spectrum.  ULrych  and  Bishop 
discuss  these  problems  and  present  an  approximate  method  to  estimate  the  filter  length.  Accord - 
ingly,  filters  one-quarter  the  length  of  the  data  were  used.  The  prediction  error  vs  filter  length 
curves  for  these  data  were  computed,  and  indicated  that  for  this  length  the  spectra  should  not  con- 
tain any  spurious  peaks  but  that  some  peaks  may  not  be  resolved.  This  situation,  rather  than  the 
reverse,  is  judged  preferable  for  this  type  of  study  since  it  is  desired  to  correlate  spectral  peaks 
that  do  not  dominate  the  spectra. 

Figure  V-7  shows  the  results  for  the  shallow  Pacific  data.  The  total  1964-74  sequence  was 
divided  into  three  equal  nonoverlapping  time  segments,  and  the  Burg  prediction  error  filter  and 
the  corresponding  spectrum  were  computed  for  each  section.  The  spectrum  for  each  segment 
shows  a peak  near  235  days.  Other  peaks,  even  though  containing  considerable  energy,  are  only 
characteristic  of  individual  segments  and  so  imply  a rather  large  degree  of  nonstationarity  in 
time  in  the  complete  sequence.  The  peak  at  235  days  contains  approximately  30  percent  of  the 
spectral  energy,  determined  by  integrating  the  spectra. 

The  largest  peak  in  the  spectrum  of  solid  earth  tides  in  the  band  considered  here  occurs  at 
10 

6 months.  Peaks  at  6 months  and  1 year  dominate  the  spectrum  of  the  rate  of  change  of  the 
1 1 

earth’s  rotation.  To  examine  if  these  periodicities  were  present  in  the  occurrence  rate  of 
earthquakes,  spectra  were  computed  on  the  whole  sequences  using  a corresponding  longer  filter 
length  (and  consequently  higher  resolving  power)  so  that  if  peaks  at  182  days  and  365  days  exist, 
then  they  would  be  resolved  from  the  peak  at  2 35  days.  Figure  V-8  gives  the  spectra  for  the 
Pacific  and  Atlantic  sequences.  As  well  as  the  235-day  peak,  a peak  at  6 months  is  present  in 
both  sequences.  The  results  for  all  ten  shallow  sequences  show  that  the  6-month  peak  is  the 
only  peak  observed  in  every  spectra.  The  mean  value  is  181.8  days,  with  a standard  deviation 
of  5.4  days.  This  peak  on  average  contains  about  5 percent  of  the  total  energy  of  the  spectra.  No 
consistent  yearly  peaks  are  observed.  The  peak  near  235  days  occurs  for  all  sequences  for  re- 
gions bordering  the  Pacific  Plate  and  the  Atlantic  Ridge.  The  mean  value  is  231.3  days,  and  the 
standard  deviation  is  18.1  days.  These  peaks  contain  an  average  of  15  percent  of  the  total  energy 
and  are  always  one  of  the  four  largest  peaks  in  the  spectrum  in  which  they  occur.  Figure  V-9 
shows  the  spectra  in  the  two  zones  outside  the  Pacific  and  Atlantic  areas.  The  6-month  peak  is 
both  time  and  space  stationary  for  these  regions;  it  accounts  for  approximately  20  percent  of  the 
energy  in  the  spectrum  of  the  occurrence  rate  in  the  Western  Mediterranean  to  China  seismic 
zone,  and  for  about  5 percent  of  the  energy  for  the  shallow  Pacific  and  Atlantic  sequences. 

In  summary,  in  the  band  2 months  to  2 years,  approximately  one -third  of  the  energy  in  the 
earthquake  occurrence  rate  spectrum  is  time  and  space  stationary,  nonrandom,  and  periodic. 

For  continental  seismicity,  the  period  is  6 months.  For  oceanic  seismicity,  the  period  is  near 
2 35  days  (approximately  7.6  months)  with  a smaller  amount  of  energy  at  a period  of  6 months. 

No  coherent  energy  at  a period  of  1 year  is  observed  for  either  the  oceanic  or  continental  seis- 
mic occurrence  rates.  , 
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Fig.  V-l.  Frequency-magnitude  curves 
for  three  tectonic  environments 


Fie.  V-2.  b-mK  plots  for  three  curves 
inFig.V-1.  b 


Fig.  V-3.  Comparison  of  frequency- 
magnitude  curves  for  Northwest  Pacific 
for  activity  maximum  of  1964-6  5 and 
period  1966-74. 
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1964  1966  1968  1970  1972  1974 


Fig.  V-4.  Earthquake  time  series  for  margins  of  Pacific  Ocean. 
PI)E  catalog,  (a)  Histogram  of  earthquake  counts  for  40-day 
intervals  (depths  100  km,  m^>5.5);  (b)  same,  smoothed  by 

sliding  200-day  window;  and  (c)  smoothed  time  series  for  deep 
events  (depth  > 300  km,  m^  > 5.2). 
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Fig.  V-5.  Time  series  for  shallow  (depths  100  km)  events  from 
regions  bordering  Pacific  Plate,  smoothed  by  20  0-day  sliding 
window.  PDE  catalog.  Magnitude  thresholds  vary  from  mu  = 
5.2  to  5.6. 
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<964  1966  <966  <9  TO  <972  <974 


Fig.  V-6.  Time  series  for  shallow  (depths  100  km)  events  from 
several  regions.  Various  catalogs  (see  text).  Thresholds  used 
are  Japan  (M  = 5.2),  Southern  California  (Mj  = 3.5),  East  Pacific 
Rise  (m^  = 4.8),  South  America  (m^  = 5.4),  Indian  Ocean  (m^  = 4.8), 
and  Atlantic  Ocean  (m^  = 5.2). 
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Fig.  V-7.  Spectra  for  three  nonoverlapping 
time  periods  of  data  set  1,  entire  shallow 
Pacific.  Earth  section  has  strong  peak  near 
235  days. 


PERIOD  (doys) 


Fig.  V-8.  Spectra  for  data  sets  1 and  3 showing  that  Atlantic 
and  Pacific  have  global  spectral  peak  near  23  5 days  individually, 
and  consistent  6 -month  peak. 
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Fig.  V-9.  Spectra  for  sequences  4 and  5.  Dominant  period  for  Mediterranean 
through  Southern  Asia  belt  is  6 months. 
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VI.  SEISMIC  DATA  MANAGEMENT  AND  COMPUTER  SYSTEMS 


A.  SEISMIC  DATA  MANAGEMENT  SYSTEM 

Work  is  continuing  within  this  system  to  develop  file  structures  and  contents  which  will 
accurately  reflect  the  status  of  seismic  data  and  information  available  from  ARPA-supported 
seismic  facilities.  Because  much  of  this  system  is  still  in  the  design  and  implementation  stages, 
the  major  effort  has  been  in  setting  up  files  and  file  structures  which  will  eventually  receive 
most  of  the  information.  The  information  already  in  existence  for  the  more  stable  components 
of  this  system  [such  as  the  large  arrays,  the  high-gain  long-period  (HGLP)  sites,  datacomputer, 
etc.]  has  been  reduced  and  structured  within  NLS  files  (the  On  Line  System  developed  by  Stan- 
ford Research  Institute).  Most  of  the  installation  reports  from  the  HGLP  sites  have  been  en- 
tered into  files  of  NLS.  These  reports  have  been  restructured  to  take  advantage  of  the  capabili- 
ties of  the  NLS.  Also  included  in  these  files  are  the  graphic  plots  for  the  seismometer  response 
curves  at  each  HGLP  site.  Generally,  most  of  the  available  information  on  the  HGLP  sites  is 
in  the  NLS.  The  one  remaining  piece  of  information  to  be  included  concerns  the  actual  data  tapes 
produced  at  each  of  the  HGLP  sites. 

The  situation  with  the  seismic  research  observatories  (SROs)  is  another  matter.  Because 
much  of  this  program  is  still  in  the  planning  and  implementation  stages,  only  information  on  the 
observatory  hardware  and  software  systems  is  now  in  the  NLS  file  system.  This  is  possible 
because  the  hardware  and  software  systems  are  nearly  identical  for  each  observatory.  However, 
because  many  of  the  actual  site  locations  are  currently  being  negotiated,  only  those  sites  which 
are  well  along  in  their  implementation  and  installation  stages  have  been  considered  at  this  point. 
In  the  master  file  containing  data  on  all  components  of  the  IWWSS,  the  individual  SROs  are  de- 
scribed in  detail.  Each  SRO  exists  as  a separate  branch  within  this  file.  Contained  in  this 
branch  is  information  on  the  site  geology,  location,  operations,  and  operator  address.  As  more 
SROs  become  operational,  each  entry  will  be  completed  for  that  particular  site. 

Other  files  contained  in  the  NLS  system  for  IWWSS  are  continually  being  updated  to  reflect 
the  progress  of  the  total  system.  Some  of  the  recent  updates  include  the  Datacomputer  Data- 
language  Manual  O/ll.  This  manual  is  being  entered  into  NLS,  and  at  this  time  is  almost  com- 
plete. Having  this  manual  located  on-line  and  available  for  immediate  reference  should  prove 
beneficial  to  users  of  the  datacomputer,  within  the  scope  of  the  IWWSS. 

R.  M.  Sheppard 


B.  DATACOMPUTER  ACCESS  EXPERIMENT 

One  of  the  major  components  in  the  IWWSS  is  the  datacomputer  at  Computer  Corporation  of 
America  (CCA).  This  data-storage  computer  will  serve  as  the  filing  system  and  storage  device 
for  the  real  seismic  data,  as  well  as  for  processed  data  such  as  bulletins,  site  status,  etc.  One 
important  question  which  we  feel  needs  investigation  concerns  the  access  time  to  a reasonably 
large  data  base  stored  in  this  datacomputer.  One  of  the  many  files  that  will  eventually  be  stored 
in  this  computer  will  be  an  event  bulletin  file.  This  file  will  probably  become  large  within  a 
short  time  after  it  is  created.  Access  to  such  a large  file  and  performing  conditional  searches 
on  this  file  could  become  a time-consuming  process.  At  this  point,  we  are  attempting  to  obtain 
some  idea  of  what  the  access  time  and  search  times  will  be  on  a file  which  may  reasonably  rep- 
resent the  bulletin  file. 
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In  an  attempt  to  produce  a file  which  will  have  some  resemblance  to  the  bulletin  file,  we 
have  organized  an  epicenter  bulletin  list  from  both  the  USGS  epicenter  list  and  the  daily  bulle- 
tins of  both  NORSAR  and  LASA.  Combining  these  bulletin  lists  for  the  period  1900  to  mid-1974 
has  produced  a bulletin  file  of  over  113,  000  events.  With  the  help  of  CCA  personnel,  the  bulle- 
tin file  has  been  stored  in  the  datacomputer.  Access  to  this  bulletin  file  is  obtained  through  the 
program  "SMART,"  written  by  personnel  at  CCA  and  runable  at  most  TENEX  sites.  This  pro- 
gram is  being  used  as  a means  of  checking  on  access  time  into  this  large  epicenter  bulletin 
data  file. 

The  initial  use  of  SMART  was  to  check  on  the  accuracy  of  the  stored  data.  Simple  data 
checks  were  performed  to  see  if  the  data  for  various  parameters  fell  within  logical  ranges. 

For  example,  SMART  was  used  to  check  the  event  depth  for  values  <0  or  >850.  In  this  manner, 
each  event  parameter  field  was  checked  to  determine  the  initial  success  in  accurately  storing 
the  seismic  data.  In  the  complete  file  of  over  113,000  events,  only  3 events  were  found  to  be 
erroneous. 

During  the  initial  check  on  access  time  with  SMART,  the  performance  of  a simple  check  on 
a single  field  took  upwards  of  19  min.  elapsed  time  and  30  sec  of  CPU  time.  These  figures 
varied  considerably,  from  11  min.  upwards  to  19  min.  for  access  and  search  time.  Some  weeks 
later,  the  same  data  search  took  from  2 to  4 min.  Since  the  program  SMART  is  being  run  at 
Office -1  on  the  ARPANET  and  is  checking  on  data  stored  at  the  datacomputer  at  CCA,  many 
factors  must  be  considered  in  the  analysis  of  the  datacomputer  access  time  across  the  ARPANET. 
More  detailed  checks  are  now  being  performed  on  this  epicenter  file  and  the  across  time  re- 
quired to  use  such  a large  data  base.  ^ ^ Sheppard 


C.  THE  ASG  ARPANET  CONNECTION 

The  Applied  Seismology  Group  (ASG)  uses  the  Lincoln  Laboratory  IBM  370/168  as  its  main 
computer  resource  for  heavy  number  crunching  and  interactive  data  processing.  Normal  access 
to  the  system  would  be  via  IBM  2741  terminals.  These  are  slow  (15  cps)  and  cannot  handle 
graphics.  To  improve  the  situation,  the  ASG  has  bought  two  Tektronix  4014-1  Graphics  Display 
Terminals.  These  are  run  at  4800  baud  (480  cps).  Rather  than  use  leased  lines,  the  connection 
to  the  Lincoln  Laboratory  machine  is  via  the  ARPANET,  a network  of  heterogeneous  computers. 
In  addition  to  giving  us  the  desired  high-speed  connections  to  Lincoln,  it  lets  us  access  the  other 
sites  with  ease.  We  are  already  using  it  to  access  the  NLS  at  Office-1  where  all  the  information 
about  the  IWWSS  is  maintained,  and  the  datacomputer  at  CCA  which  will  at  some  point  contain 
a large  seismic  base. 

All  our  terminals  (two  4014's,  two  TI  silent  700’s,  and  a Datapoint)  are  connected  to  a 
PDP-ll/50  at  our  Carleton  Street  site.  This, 'in  turn,  is  connected  via  a special  interface  to 
the  network.  On  the  11/50  we  run  the  ELF  software  system  to  support  the  various  network  pro- 
tocols. In  the  current  configuration,  it  can  support  up  to  6 terminals  simultaneously  (at  present 
we  have  only  5 terminals).  Reliability  is  tolerable  for  an  interim  system:  3 crashes/day  on 
4 days  out  of  5;  10  crashes  on  the  fifth  day. 

In  addition  to  providing  terminal  access  to  the  ARPANET,  the  system  allows  the  users  to 
get  listings  from  the  remote  site  and  print  them  locally.  This  facility  is  used  heavily  from 
various  TENEX  sites  (Office-1,  BBN,  etc.). 

We  also  have  a PDP-ll/40  with  64K  core.  Hardware  reconfiguration  enabling  the  11/40  to 
run  ELF  is  under  way.  When  completed,  we  will  be  able  to  run  more  users  and  have  room  to 
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insert  software  to  handle  specific  ASG  problems.  The  11  /50  will  then  be  free  for  program 
development  and  direct  support  of  seismic  data  processing. 

Y.  Peduel 


D.  PDP-7  CONSOLE  EXPANSION 

As  our  research  effort  emphasis  goes  from  a few  large  arrays  to  networks  of  stations, 
certain  limitations  of  the  Seismic  Data  Analysis  Console"!  and  of  the  data  formats  which  it  can 
use  have  become  apparent.  We  have  therefore  undertaken  modifications  and  software  develop- 
ment which  will  allow  us  to  use  the  Console  system  more  effectively  with  seismograms  from  a 
network  of  stations.  This  will  allow  us  to  utilize  the  power  and  features  of  the  Console  for  net- 
work data  without  developing  a completely  new  system  for  the  purpose. 

The  planned  changes  to  the  Console  are: 

(1)  Modify  the  CRT  display  to  show  more  identification  information  for  indi- 
vidual traces.  This  has  been  completed. 

(2)  Allow  the  user  to  operate  a knob  to  time-shift  a single  seismic  trace  on 
the  display  without  changing  other  displayed  traces,  and  to  selectively 
change  the  gain  of  individual  traces.  The  required  software  changes  to 
achieve  this  have  been  identified,  and  coding  is  in  progress. 

(3)  Provide  for  efficient  event-oriented  initialization  of  the  Console  from  a 
new  format  data  tape  containing  network  seismograms.  In  addition  to 
the  seismogram  itself,  the  new  tape  format  will  contain  all  station  and 
instrument  information  required  by  the  initialization  program  for  each 
seismogram.  The  structure  of  the  Console  system  will  allow  us  to  write 
the  new  initialization  capability  in  Fortran  once  a small  number  of  as- 
sembly language  programs  are  written.  Neither  the  design  nor  assembly 
language  subroutines  are  currently  completed. 

(4)  Provide  for  dumping  the  data  of  a Console  session,  including  station  iden- 
tification, instrument,  and  time  information  for  each  individual  channel, 
on  tape  for  future  processing.  A simple  program  to  do  this  has  been 
written.  It  generates  a tape  which  subsequently  can  be  used  to  re-initialize 
the  Console  with  a single  command.  With  some  minor  modifications  this 

is  sufficient  for  the  immediate  future,  but  eventually  the  capability  to  out- 
put in  a more  standard  form  for  use  by  other  programs  will  be  required. 

The  new  format  input  tapes  to  the  Console  mentioned  above  will  be  quite  general  and  not  just 
limited  to  Console  application.  They  will  contain  four  distinct  sections,  namely:  (1)  Type, 

(2)  Description,  (3)  Title,  and  (4)  Body.  The  four  of  these  together  will  constitute  a Data  Unit 
(DU).  The  Body  of  a DU  will  contain  actual  data.  What  is  stored,  how  it  is  stored,  and  how  it 
can  be  manipulated  are  determined  by  the  Type  and  Description  sections.  The  Description 
section  will  specify  the  logical  structure  of  the  data  items  in  the  Body,  while  the  Type  section 
will  specify  the  storage  device  and  the  physical  representation  of  the  items.  The  Title  will  be 


t P.  L.  Fleck  and  L.  J.  Turek,  "A  Seismic  Data  Analysis  Console,1’  Technical  Report  495,  Lin- 
coln Laboratory,  M.l.T.  (18  January  1972),  DDC  AD-740604. 
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a short  section  containing  creation  date  and  some  alphanumeric  text  supplied  by  the  user  to  help 
him  identify  the  data  in  the  DU.  A program  using  a DU  must  first  check  the  Type  and  Descrip- 
tion to  determine  whether  it  can  handle  that  DU.  The  representation  of  the  Type  and  Description 
sections  is  a function  of  the  storage  device  containing  the  DU.  Thus,  when  a DU  is  accessed  we 
know  automatically  how  to  read  these  sections  because  we  know  what  device  we  are  using. 

Our  first  planned  use  of  Data  Units  is  to  input  multi-site  waveform  data  to  the  Console.  In 
particular,  the  SROs  and  I1GLP  stations  are  of  interest.  Such  data  will  be  input  to  the  Console 
via  a single  family  of  DUs  which  have  a common  Description.  The  Description  will  specify  that 
the  Body  is  a stream  of  objects,  each  object  consisting  of  a header  giving  start  times,  station 
name,  etc.,  and  an  array  of  waveform  samples.  Although  not  specified  by  the  Description,  the 
objects  will  be  ordered  according  to  the  start  times.  Long  waveforms,  say  more  than  a thousand 
samples,  will  simply  be  broken  into  a number  of  objects.  The  Console  initialization  program 
will  work  from  a condition  table  listing  the  station,  instrument,  and  time  intervals  of  interest, 
and  then  scan  through  the  DU  transferring  to  the  Console  all  data  that  satisfy  the  given  conditions. 
This  approach  gives  great  flexibility  in  the  way  the  user  can  request  data,  since  the  condition 
table  can  be  set  up  in  a variety  of  ways.  For  example,  suppose  an  initialization  program  to  get 
short-period  P-wave  data  based  upon  predicted  P-arrival  times  for  an  event  is  required  but  docs 
not  exist.  It  will  be  necessary  only  to  write  a Fortran  program  to  create  a condition  table  con- 
taining the  computed  arrival  times  at  the  stations  of  interest. 

A program  is  being  written  to  convert  a single  9-track  SRO  tape  into  a DU  on  a 7-track  tape 
which  can  be  read  by  the  PDP-7  computers  which  support  the  Console.  In  addition,  a program 
to  merge  DUs  from  different  stations  on  the  basis  of  signal  start  times  is  being  written.  This 
merge  is  required  to  obtain  a library  of  DU  tapes  with  all  network  data  for  a time  period  on  a 
single  tape.  This  can  be  done  routinely  and  means  that  subsequent  Console  users  generally  will 
need  to  mount  only  one  tape  to  get  all  available  data  for  an  event. 

R.  T.  Lacoss 

L.  J.  Turck 

M.  F.  O’Brien 


E.  SOFTWARE  LIBRARIES 

Seismic  research  activities  are  aided  considerably  by  general  and  application  program 
libraries  which  are  adequately  documented  and  easy  to  use.  In  the  past  we  have  developed,  and 
we  continue  to  expand,  a very  extensive  library  for  our  PDP-7  computers.  This  library  is  on- 
line and  always  available  to  users  and  programs.  In  addition,  a smaller  library  has  been  main- 
tained for  batch  processing  on  the  Lincoln  IBM  370/168.  This  has  been  limited  largely  to  rou- 
tines to  read  data  tapes  and  to  perform  functions  such  as  high-resolution  spectral  analysis  which 
require  considerable  computing  power. 

Our  ARPANET  connection  to  the  Lincoln  370/168  has  enhanced  our  access  to  that  facility 
and  its  utility  to  us.  To  make  the  best  use  of  this,  we  have  undertaken  to  significantly  expand 
and  reorganize  our  program  libraries  on  that  system.  Although  some  of  this  will  involve  changes 
in  the  library  available  for  batch  processing,  the  major  effort  is  directed  toward  a new  library 
for  interactive  use. 

Work  on  the  new  library  for  time-sharing  users  of  the  IBM  370/168  has  just  been  started. 

The  first  stage  is  primarily  the  conversion  of  existing  programs  and  generation  of  proper  docu- 
mentation for  their  use  in  the  new  environment.  Since  many  programs  actually  are  interfaces 
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to  data  bases  such  as  names  and  locations  of  seismic  stations,  travel -time  curves,  and  digitized 
maps,  the  new  library  will  contain  such  data  as  well  as  actual  programs.  A routine  to  obtain 
seismic  station  locations  from  their  names  is  now  available.  Other  basic  routines  currently 
available  to  users  on  our  PDP-7  systems  will  be  added  shortly. 

R.  M.  Sheppard 
L.  Lande 
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